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COMPUTATION  OF  PROPAGATION  SPEED  AND  REFLECTION  OF 
AXIALLY  SYMMETRIC  WAVES  IN  COMPOSITE  CYLINDERS ,  WITH 
APPLICATION  TO  IMPEDANCE  TUBE  AND  CALIBRATOR 


1.0.  INTRODUCTION 

Various  instruments  used  at  the  Underwater  Sound  Reference  Detachment  of 
the  Naval  Research  Laboratory  (NRL-USRD)  are  based  on  acoustic  waveguide 
technology.  Examples  are  the  impedance  tube  [1],  inertial  calibrators  [2,3], 
and  the  Low-Frequency  Facility's  systems  [4,5].  A  common  feature  of  these 
instruments  is  the  creation  of  an  acoustic  field  with  well-defined  properties 
in  a  limited  space,  as  contrasted  with  free-field  measurements.  The  field 
description  of  this  class  of  instruments  is  based  on  waveguide  theory.  In 
this  report  an  analysis  of  the  waveguide  aspects  of  these  instruments  is 
presented.  Discussion  and  documentation  of  the  computer  programs  used  in  the 
data  reduction  are  given.  The  correspondence  between  predictions  from  the 
analysis  and  the  data  will  be  a  measure  of  the  extent  to  which  simple 
waveguide  theory  is  adequate  in  describing  the  phenomena. 

Waveguide  theory  is  well  established  and  described  1 6 , 7  J .  Some  of  its 

features  are  given  in  this  report  to  serve  as  a  basis  for  the  understanding  of 
the  computer  programs  and  quantitative  conclusions.  It  is  assumed  that  the 
wave  propagation  takes  place  in  one  dimension  only.  The  field  perpendicular 
to  this  direction  is  spatially  limited  by  shape  and  properties  of  physical 
bodies.  This  does  not  exclude  the  possibility  that  particle  velocity  and 
pressure  may  be  functions  of  more  than  one  coordinate.  It  will  be  assumed, 
though,  that  at  interfaces  perpendicular  to  the  direction  of  propagation,  the 
boundary  conditions  are  only  imposed  on  averaged  axial  velocity  and  averaged 
pressure.  Usually,  one  applies  waveguide  theory  to  cases  where  the  cross- 
sectional  dimension  is  small  compared  with  the  wavelength;  a  consequence  is 
that  the  radial  variation  of  the  field  variables  is  small.  The  analysis 
developed  in  this  report  does  not  account  for  the  influence  of  finite  cylinder 
length  and  its  attending  termination  on  the  wave  pattern.  This  influence  will 
be  less  the  larger  the  ratio  of  length  to  diameter  of  the  various  cylindrical 
elements  of  the  system. 

The  hierarchy  of  computer  programs  documented  here  may  be  understood  by 
reference  to  the  sketches  of  impedance  tube,  Fig.  1,  and  free  surface 
calibrator.  Fig.  2.  For  the  cylindrical  sample  in  the  impedance  tube,  Fig.  1, 
one  wants  to  evaluate  the  propagation  speed  in  the  axial  direction  given  the 
values  of  a  pair  of  elastic  moduli  of  the  solid  material.  In  general, 
propagation  speed  and  elastic  moduli  are  complex,  accounting  for  the  existence 
of  dissipation  of  acoustic  energy  in  the  material.  The  reflection  of  a  plane 
wave  from  the  front  face  of  the  sample  may  then  be  computed.  The  propagation 
speed  in  the  water  will  be  influenced  by  the  lark  of  rigidity  of  the  cylinder 
wall.  This  problem  is  more  serious  in  the  case  of  the  calibrator  where  the 
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wall  thickness  is  considerably  less  than  that  of  the  impedance  tube.  This 
influence  is  numerically  determined.  The  sequence  of  computations  in  data 
analysis  of  the  impedance  tube  is  just  the  opposite  of  the  one  sketched.  One 
experimentally  determines  the  complex  reflection  coefficient.  From  this  a 
complex  propagation  speed  is  inferred  and,  in  turn,  this  propagation  speed  is 
related  to  the  elastic  moduli  of  the  solid  material.  This  inverse  procedure 
may  be  effected  by  a  combination  of  root-search  techniques  and  the  programs 
presented  here. 


Fig.  1  -  Section  of  acoustic  impedance  tube. 


Fig.  2  -  Principle  of  free  surface  acoustic  calibration  device. 

In  this  report  three  computer  programs  are  described  and  documented. 
Program  RTUBE  computes  the  complex  propagation  speed  in  the  sample  in  the 
impedance  tube  from  the  measured  complex  reflection  coefficient.  Program 
WGUIDE,  with  various  options,  determines  real  propagation  speed  in  the  sample 
in  the  impedance  tube  and  in  the  water  in  the  calibrator  from  a  pair  of  real 
elastic  moduli  characterizing  the  material  of  the  sample  and  tube  wall, 
respectively.  Program  IMPED  computes  the  complex  wave  speed  from  two  complex 
elastic  moduli  for  the  sample  in  the  impedance  tube. 


In  Che  next  section  Che  theory  for  reflection  of  the  signal  in  the 
impedance  tube  is  given.  Then  the  theory  of  composite  waveguides  is  developed 
and  the  various  dispersion  relations  are  derived  applicable  to  the  boundary 
conditions  for  impedance  tube  and  calibrator.  Examples  are  given  for  wave 
propagation  in  samples  of  steel  and  rubber  in  the  impedance  tube,  and  for  the 
influence  of  the  wall  of  the  calibrator.  The  detailed  discussion  and 
documentation  of  the  computer  programs  is  given  in  Appendix  A.  Appendix  B 
deals  with  the  application  bf  a  correction  .factor  to  the  measurement  in  the 
impedance  tube  that  accounts  for  the  lack  of  perfection  of  the  reflector. 

An  important  conclusion  of  this  study  is  the  observation  that  the  sound 
speed  in  the  sample  in  an  impedance  tube,  ignoring  end  effects,  is  not 
necessarily  the  dilatational  wave  speed  but  a  function  of  the  size  of  the 
fluid-filled  gap  between  sample  and  tube  wall,  and  the  elastic  properties  of 
the  sample  material.  This  effect  may  be  utilized  to  obtain  more  information 
◦n  the  elastic  moduli  by  varying  the  sample  radius,  and  also  by  not  admitting 
fluid  to  fill  this  gap. 

2.0.  REFLECTION  COEFFICIENT  IN  IMPEDANCE  TUBE 

It  is  assumed  that  a  wave  propagates  in  the  axial  direction  identified 
with  the  z  coordinate  and  that  the  wave  pattern  is  periodic  in  z.  Thus  the 
field  variables  will  contain  a  factor  exp  i(wt  -  kz),  where  w  is  the  angular 
frequency  and  k  the  wave  number;  the  phase  speed  c  is  related  by  u  ■  kc. 
Because  of  the  geometry,  it  is  assumed  that  the  field  variables  do  not  depend 
on  the  azimuth  angle.  The  pressure  and  axial  particle  velocity  are,  in 
general,  functions  of  the  radial  distance  r.  For  the  sake  of  the  computation 
of  the  reflection  coefficient,  it  is  assumed  that  one  may  represent  the  wave 
as  a  plane  wave  with  average  values  for  pressure  and  velocity.  The  sample  of 
length  d  is  backed  by  a  reflector  with  acoustic  impedance  Zr.  One  can  show 
[8]  that  the  input  impedance  at  the  front  end  of  the  sample  is 


Z^  +  i(pgc/S)  tan  kd 
(pgc/S)  ^  C/S)  ♦  i  Z  tan  kd  ’ 


(1) 


where  p  is  the  density  of  the  sample 


c  is  the  sound  speed  in  the  sample 
S  is  the  cross-sectional  area 
k  is  the  wave  number  in  the  sample. 


If  one  applies  this  same  formula  to  the  reflector  behind  the  sample, 
where  the  impedance  of  the  air  backing  the  reflector  may  be  ignored,  one  sees 
that  the  impedance  of  the  reflector,  Z  ,  is  given  by  Zf  *  i(pc)p  tan  kr4, 
where  the  subscript  r  refers  to  the  reflector,  and  4  is  the  length  of  the 
reflector.  In  practice,  one  chooses  4  such  that  k4  ■  n/2  for  some 
intermediate  frequency  and  then  the  impedance  Zf  is  virtually  infinite  for  a 
reasonably  large  bandwidth. 


The  reflection  coefficient  (complex)  r  of  the  sample  is  defined  as  the 
ratio  of  the  amplitude  of  the  reflected  to  that  of  the  incoming  wave. 
Equating  impedances  at  the  interface  leads  to  the  relation  between  r  and  the 
wave  number  k  according  to 


1  ♦  r 
1  -  r 


-  -i ( p  / p  )  (k  d/kd)  cot(kd) 
so  o 


where  pQ  is  Che  density  of  and  kQ  Che  wave  number  in  Che  fluid.  The  program 
RTUBE  determines  the  complex  root  kd  of  this  equation  from  the  measured 
complex  value  of  r. 

The  data  reduction  incorporated  in  program  RTUBE  is  based  on  Eq.  (2) 
which  is  derived  under  the  assumption  that  the  reflector  is  ideal;  i.e.,  that 
the  amplitude  and  phase  of  the  reflected  pressure  wave  are  the  same  as  those 
of  the  incident  wave.  To  account  for  the  lack  of  rigidity  of  the  reflector, 
and  to  avoid  measurement  of  absolute  magnitude  and  phase,  the  wave  reflected 
in  the  tube  with  the  sample  removed  is  measured.  The  reflection  coefficient  r 
is  computed  according  to  r  ■  (r  /r  )exp(-2ik  d),  where  r„_  is  the 
reflection  with  the  sample  in  place,  and  rgt  is  the  reflection  without  the 
sample.  The  exponential  factor  accounts  for  the  fact  that  the  sample  is 
replaced  by  a  column  of  water  with  a  given  phase  shift.  This  r  is  only 
approximately  correct;  an  additional  correction  factor  is  described  in 
Appendix  B. 

3.0.  THEORY  Of  HAVE  PROPAGATION  IN  INFINITE  CYLINDERS  AND  CYLINDRICAL  SNELLS 

In  an  elastic  solid  of  infinite  extent  there  are  two  different  types  of 
wave  motion.  One  is  an  irrotational  wave,  usually  denoted  as  a  dilatational 
wave,  and  the  other  is  an  incompressible  motion,  or  shear  wave. 

Mathematically,  they  can  be  described  by  a  scalar  and  vector  potential  A  and 
respectively.  In  the  axially  symmetric  situation  of  waveguides,  only  the 
azimuthal  component  of  the  vector  potential  is  needed,  henceforth  indicated  by 
the  function  H(r,z). 

The  general  relations  used  below  can  be  found  in  Ref.  9.  The  particle 
displacement  u  is  given  in  terms  of  the  potentials  by 
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pertinent  to  the  analysis  for  the  case  of  axial  symmetry  are 
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and  X  and  y  are  Lam£  constants. 


The  equations  of  motion  are 


3t  3t 
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3r  +  3z  +  r  Trz  Ps  3t2 


(5) 


By  combining  Eqs.  (3)  through  (5)  one  derives  the  following  wave  equations  for 
the  potentials  <p  and  H 

c  2A$  *  — —  (6) 

d  3t2 

and 

2 

c  2  (AH  -  H/r2)  =  .  (7) 

8  3t2 

where  c^  is  the  wave  speed  for  dilatational  waves,  c^2  =  (X  +  2y)/ps, 
c  is  the  wave  speed  for  shear  waves,  and  c  2  =  u/p  ;  A  is  the  Laplacean 


operator  for  cylindrical  coordinates  A 
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If  one  assumes  solutions  for  $  and  H  in  the  form  of  a  wave  traveling  in 
the  z  direction,  according  to 


and 


,,  ,  i(ut  -  kz) 

<(>  ■  f(r)  e 

„  ,  .  i(ut  -  kz)  , 

H  -  g(r)  e 


(8) 


the  functions  f(r)  and  g(r)  are  solutions  of 


f"  +  f'/r  -  q2f  -  0 


(9) 


and 


g"  +  g'/r  -  (s2  +  Vr2)g  *  0  ,  (10) 

where  q2  «*  k2  -  k2  and  s2  =  k2-k2,k,  and  k  are  the  wave  numbers  for 
dilatational  and  shear  waves.  In  terms  or  the  viriables  q'  and  s',  defined  by 
(q')2  =  -q2  and  (s')2  *  -s2,  the  general  solutions  to  Eqs.  (9)  and  (10)  are 
f  *  A  J  (q*r)  +  B  Y  (q'r)  and  g  =■  C  Ji(s'r)  +  D  Yi(s'r).  J  is 
the  Bessel  function  of  the  first  kind  with  order  n,  and  Yn  isnthe  Bessel 
function  of  the  second  kind  with  order  n. 

The  compress ional  wave  in  the  fluid  is  represented  by  a  velocity 
potential  which  is  a  solution  of  the  wave  equation 

c  2A$  -  - —  (H) 

°  °  3t2 

where  cQ  is  the  wave  speed  in  the  fluid. 

The  velocity  components  U  and  U  of  the  particle  velocity  in  the  fluid 
3*o  3*o 

follow  from  Ur  *  —  and  Uz  *  .  The  solution  of  Eq.  (11)  representing  a 

traveling  wave  in  the  z  direction  is 


*  «  IE  J  (q'r)  +  F  Y  (q'r)]  ei(u>t  "  kz) , 
o  o  o  o  o 

where  (q ' ) 2  ■  k  2  -  k2,  k  is  the  wave  number  of  the  compressional  wave  in  the 
fluid.  Application  of  the  boundary  conditions  leads  to  the  dispersion 
relations  for  the  phase  speed  c.  For  a  wave  in  an  infinite  cylinder  in 


vacuum,  the  analysis  has  been  performed  by  Pochhammer  and  Chree  [10].  A 
comprehensive  set  of  calculations  of  the  phase  speed  based  on  the  Pochhamraer- 
Chree  solution  were  performed  by  Bancroft  [11]. 

4.0.  PROPAGATION  SPEED  OF  AXIALLY  SYMMETRIC  HAVES  IN  COMPOSITE  CYLINDERS 

The  determination  of  phase  speed  in  composite  cylinders  is  accomplished 
by  imposing  the  pertinent  boundary  conditions  on  the  stresses  and  particle 
velocities  derived  in  Section  3.  An  overview  of  the  various  optior  in  the 
computer  program  is  presented  in  Fig.  3. 

The  physical  situation  in  the  impedance  tube  is  represented  it 
Option  1.  A  solid  cylinder  is  surrounded  by  a  thin  layer  of  fluid  d  the 
center  wall  of  the  tube  is  assumed  to  be  rigid.  For  comparison,  a  outer 

boundary  is  considered  in  Option  2,  although  this  case  can  hardly  i  alized 
in  practice.  The  calibrator  situation  is  represented  by  Option  3,  where  the 
influence  of  a  nonrigid  wall  on  the  fluid  wave  speed  is  evaluated.  Options  4 
and  3  reflect  the  situation  where  a  soft  coating  is  applied  to  the  inside  of 
the  calibration  tube  in  order  to  slow  down  the  wave  in  the  liquid;  in 
Option  4  the  coating  is  cemented  to  the  wall,  in  Option  5  it  is  free  to  move 
tangentially  at  the  wall.  In  Option  6  the  wave  speed  in  a  fluid  cylinder  with 
rigid  boundary  is  computed,  in  Option  7  the  wave  speed  in  a  fluid  cylinder 
with  free  boundary  is  computed. 


OPTION  3  OPTION  * 


Fig.  3  -  Cross  section  of  waveguides  in  situations  corresponding 
to  options  in  program  WGUIDE. 
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In  the  real  variable  program  WGUIDE,  the  variables  q' ,  s',  and  q'  may  be 
either  real  or  imaginary.  In  order  to  avoid  the  necessity  of  complex0 
calculations,  the  transformations  q'  ♦  iq;  s'  ♦  is,  and  q'  iq  are 
applied  whenever  imaginary  values  appear.  The  Bessel  functions  are 
transformed  according  to  the  identities  J0(i*>  ■  IQ(x),  J^(ix)  ■  il^(x), 

YQ(ix)  *  -K  (x),  and  Yj(ix)  -  iKj(x).  IQ(x)  K^x)  are  “°dified  Bessel 
functions  [12 J. 


4.1.  Have  Speed  in  Sample  in  Impedance  Tube 


The  following  set  of  boundary  conditions  is  applied  to  the  field 
variables  in  the  solid  cylinder  and  surrounding  fluid  if  the  outer  wall  is 
rigid  (Option  1). 

a.  The  normal  stress  Trr  in  the  solid  is  equal  to  the  negative  of 
the  pressure  p  in  the  fluid  at  the  surface  r  ■  a,  where  a  is  the 
radius  of  the  cylinder.  The  pressure  follows  from 


it 

3r 


Thus  p 


-P 


3<t> 

c 

o  3t 


at  r 


a. 


b.  The  tangential  stress  Tfz  is  zero  in  the  solid  for  r  *  a. 

c.  The  normal  component  of  the  particle  velocity  is  continuous  at 
the  interface  between  solid  and  fluid  at  r  ■  a. 

d.  The  normal  component  of  the  particle  velocity  in  the  fluid  Ur  is 
zero  at  r  *  b. 

The  coefficients  of  the  amplitudes  A,  C,  E,  and  F  may  be  shown  in  a 
matrix,  Table  1.  Equating  the  determinant  value  of  this  matrix  to  zero  yields 
the  dispersion  relation  for  the  wave  speed. 

If  the  outer  boundary  is  assumed  free  (Option  2\  the  fourth  boundary 
condition  is  replaced  by  one  where  the  pressure  in  the  fluid  is  zero  for 
r  *  b.  The  resulting  matrix  of  the  coefficients  is  shown  in  Table  2. 
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4.2.  Have  Speed  in  Calibrator 

Option  3 


For  the  case  of  the  calibrator,  one  assumes  that  the  fluid  cylinder  has  a 
radius  a  and  is  bounded  by  a  solid,  but  not  rigid,  wall  with  outer  radius  b. 
The  following  set  of  boundary  conditions  applies  to  the  field  variables  in 
fluid  and  cylinder  wall. 

a.  The  negative  value  of  the  pressure  p  is  equal  to  the  normal 
stress  in  the  solid  at  r  ■  a. 

b.  The  tangential  stress  in  the  solid  is  zero  at  r  *  a. 

c.  The  normal  stress  in  the  solid  is  zero  at  r  *  b. 

d.  The  tangential  stress  in  the  solid  is  zero  at  r  =  b. 

e.  The  normal  component  of  the  particle  velocity  in  the  fluid  is 

equal  to  that  in  the  solid  at  r  *  a. 

The  resulting  coefficient  matrix  for  amplitudes  A,  B,  C,  D,  and  F  is  shown  in 
Table  3. 

It  was  desirable  to  analyze  the  situation  where  a  soft  coating  was 
applied  to  the  inside  of  the  calibrator  with  the  purpose  of  reducing  the 
effective  wave  speed  in  the  liquid  of  the  calibrator.  If  this  coating  is 
cemented  to  the  wall,  the  following  boundary  conditions  apply: 


Option  4 


hi 


a.  The  negative  value  of  the  pressure  p  is  equal  to  the  normal  stress  in 
the  solid  at  r  »  a. 

b.  The  tangential  stress  in  the  solid  is  zero  at  r  *  a. 

c.  The  normal  component  of  the  velocity  is  zero  at  r  *  b. 

d.  The  tangential  component  of  the  particle  velocity  is  zero  at  r  “  b. 

e.  The  normal  component  of  the  particle  velocity  is  continuous  at  r  ■  a. 

The  resulting  matrix  of  the  coefficients  of  these  boundary  conditions  is 
shown  in  Table  4. 


« 
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Matrix  of  coefficients  of  boundary  conditions  for  coating  cemented  to  calibrator  wall 


If  the  coating  Is  not  cemented  to  the  wall,  but  is  free  to  move  in  a 
tangential  direction,  condition  (d)  is  replaced  by  the  requirement  that  the 
tangential  stress  in  the  solid  is  zero  at  r  -  b.  The  resulting  matrix  is 
given  in  Table  5. 


Table  5  -  Matrix  of  coefficients  of  boundary  conditions  for  coating  not  cemented  to  calibrator  wall. 


Empty  cylinder 


The  phase  speed  in  an  empty  solid  cylinder  may  be  readily  derived  from 
the  matrix  in  Table  2  by  omitting  the  fifth  row  and  fifth  column  and  setting 
the  determinant  value  of  the  remaining  4*4  matrix  equal  to  zero. 

Other  phase  speeds  of  Interest  may  be  obtained.  In  a  fluid  column  with 
free  boundaries,  the  phase  speed  is  given  by 

J  (q*  a)  -  0  or  (c/c  )2  -  1/[1  -  j*/(k  a)2]  (12) 

o  o  o  o,m  o 


where  m  ■  1,2,...  and  jQ  m  are  the  zeros  of  the  JQ  Bessel  function.  Thus, 
there  is  a  cutoff  frequency  given  by 


u 

o 


Jo,lCo/a 


(13) 


below  which  no  wave  propagation  is  possible  in  a  free  fluid  cylinder. 
Similarly,  the  phase  speed  of  waves  in  a  fluid  cylinder  bounded  by  rigid  walls 
follows  from 


(V>Ji(q0*)  " 0  (14) 

or 

(c/cQ)2  -  1/[1  -  j12m/(koa)2l  (15) 

where  m  »  0,1,2,...  and  m  are  the  zeros  of  the  Bessel  function 
(including  Ji fo  “  °) • 

4.3.  Complex  Have  Speed  in  Sample  In  Impedance  Tube 

The  measurements  in  the  impedance  tube  consist  of  a  complex  reflection 
coefficient.  Thus,  it  is  important  to  determine  the  complex  wave  speed  that 
is  connected  with  a  pair  of  complex  elastic  moduli.  This  is  accomplished  by 
the  program  IMPED.  Its  dispersion  relation  is  the  same  as  in  WGUIDE  option  1, 
but  the  roots  are  found  by  a  complex  root  finder  described  in  Ref.  13. 


5.0.  RESULTS  AMD  DISCUSSION 
5.1.  Impedance  Tube 

The  original  Incentive  to  develop  the  WGUIDE  program  was  the  need  to 
determine  exactly  which  propagation  speed  is  measured  in  the  impedance  tube. 

It  appears  to  be  tacitly  assumed  that  the  dilatational  speed  is  the  quantity 
measured.  For  instance,  in  Ref.  14  the  sound  speed  is  graphed  without  further 
qualification;  from  the  text  one  may  infer  that  it  is  the  dilatational  wave 
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speed.  The  present  analysis  shows  that  it  depends  on  the  circumstances  as  to 
what  speed  applies  for  the  sample  in  the  impedance  tube.  This  is  shown  in  a 
sequence  of  figures.  The  first  set  shows  the  dispersion  relation  in  a  sample 
of  steel  with  a  density  of  7700  kg/m3,  Young's  modulus  19.5xl010  Pa, 
shear  modulus  8.3x10*®  Pa,  Poisson's  ratio  0.37  (data  from  Ref.  15).  The  tube 

is  filled  with  water  with  density  998  kg/m3  and  wave  speed  1481  m/s.  The 

abscissa  kga  is  a  dimensionless  measure  of  the  frequency,  with  kg  the  wave 
number  for  the  dilatational  wave  in  the  solid  and  a  the  radius  of  the 
cylindrical  sample;  the  ordinate  is  the  ratio  of  the  wave  speed  to  the 
dilatational  wave  speed.  For  large  values  of  kga  the  wave  speed  approaches 
the  Rayleigh  wave  speed.  Figure  4  shows  that  for  a  sample  in  vacuum  at  low 
frequency  the  bar  speed  is  approached,  as  expected.  With  water  in  the  space 
between  the  tube  wall  and  the  sample,  radii  a  and  b,  respectively,  the  wave 
speed  in  the  sample  at  low  frequency  increases  with  decreasing  value  of  b/a  as 
seen  in  Figs.  5  and  6  and  is  indistinguishable  from  the  dilatational  wave 
speed  (c/cd  ■  1)  for  b/a  *  1.001,  Fig.  7.  Notice  that  there  is  a  "fluid  mode" 
with  low  propagation  speed  in  addition  to  various  modes  determined  by  the 
sample.  The  same  behavior  is  shown  by  hard  rubber  as  sample  material,  density 
1100  kg/m3,  shear  modulus  0.1x10*®  Pa,  bulk  modulus  0.5x10*®  Pa.  Here, 
however,  a  thicker  layer  of  fluid  is  able  to  constrain  the  material  to 
manifest  dilatational  wave  speed  than  for  steel  (see  Figs.  8  through  10). 

Even  for  b/a  ■  1.2,  the  speed  is  close  to  the  dilatational  wave  speed  (Fig.  9) 

The  value  for  b/a  is  about  1.03  for  the  NRL-USRD  impedance  tube. 


Fig.  4  -  Relative  wave  speed  c/cg  as  a  function  of  dimensionless  frequency 
kga  for  steel  cylinder  in  vacuum.  Vertical  arrows  indicate 
asymptotes;  Cg  is  the  dilatational  speed,  cY  is  bar  speed. 
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RELATIVE  WAVE  SPEED  c/c, 


SAMPLE  MODES 


FLUID  MOOE 


DIMENSIONLESS  FREOUENCY  kdo 

Relative  wave  speed  c/c.  as  a  function  of  dimensionless  frequency 
k^a  for  steel  cylinder  in  water,  bounded  by  a  concentric  rigid 
wall.  Ratio  of  the  cylinder  radii  is  1.05.  Vertical  arrows  indi¬ 
cate  asymptotes;  c^  is  the  dilatational  speed,  c„  is  the  bar  speed 


DIMENSIONLESS  FREQUENCY  kda 


z 


5.0 


Relative  wave  speed  c/c^  as  a  function  of  dimensionless  frequency 
k^a  for  a  steel  cylinder  in  water  bounded  by  a  concentric  rigid 
wall.  Ratio  of  cylinder  radii  is  1.001.  Vertical  arrows  indicate 
asymptotes;  d^  is  the  dilatational  wave  speed,  Cy  is  the  bar  speed. 


DIMENSIONLESS  FREQUENCY  kja 


Relative  wave  speed  c/c^  as  a  function  of  dimensionless  frequency 
k^a  for  a  rubber  cylinder  in  vacuum.  Vertical  arrows  indicate 
asymptotes;  cd  is  the  dilatational  wave  speed,  Cy  is  the  bar  speed 


DIMENSIONLESS  FREQUENCY  Lda 


9  -  Relative  wave  speed  c/c^  as  a  function  of  dimensionless  frequency 
k^a  for  a  rubber  cylinder  in  water,  bounded  by  a  concentric  rigid 
wall.  Ratio  of  cylinder  radii  is  1.2.  Vertical  arrows  indicate 
asymptotes;  is  the  dilatational  wave  speed,  Cy  is  the  bar  speed 
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DIMENSIONLESS  FREQUENCY  kda 


Fig.  10  -  Relative  wave  speed  c/c^  as  a  function  of  dimensionless  frequency 
k^a  for  a  rubber  cylinder  in  water,  bounded  by  a  concentric  rigid 
wail.  Ratio  of  the  cylinder  radii  is  1.05.  Vertical  arrows 
indicate  asymptotes;  c^  is  the  dilatational  wave  speed,  Cy  is  the 
bar  speed. 


Generally,  for  elastomers,  the  measured  wave  speed  is  correctly 
identified  with  the  dilatational  wave  speed  since  the  compressibility  and 
density  of  most  rubbers  are  close  to  that  of  water,  but  this  is  not  true  for 
more  rigid  materials.  A  single  measurement  is  not  sufficient  to  determine  the 
•  pair  of  elastic  moduli  needed  to  describe  the  material;  by  preventing  the  fill 

"  fluid  to  reach  the  gap  between  the  sample  and  the  tube  or  by  varying  the 

diameter  of  the  sample,  one  may  measure  values  for  different  combinations  of 
the  moduli  from  which  the  desired  information  may  be  extracted. 


5.2.  Calibrator 

Two  issues  of  interest  were  investigated  with  the  WGUIDE  program.  In  the 
first  place,  one  may  consider  the  influence  of  the  walls  on  the  sound  speed  in 
the  water.  One  expects  that  the  less  rigid  the  wall  the  more  the  speed  will 
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be  reduced.  As  an  example,  the  dispersion  relation  was  computed  for  the  case 
of  the  G40  calibrator  developed  at  NRL-USRD.  The  following  data  apply: 
a  =  9.68  cm,  b  ■  10.95  cm,  the  material  of  the  wall  is  aluminum  with  a  density 
of  2700  kg/m^,  shear  modulus  2.4x10*0  Pa,  bulk  modulus  7.5x10*0  Pa 
and  a  dilatational  wave  speed  of  6300  m/s  (Ref.  15).  The  results  are  shown  in 
Fig.  11.  It  is  seen  that  at  low  frequency,  (k^a  *  0.1,  which  corresponds  to 
1035  Hz),  the  ratio  of  the  speed  in  the  tube  to  that  in  water  is  0.83.  Thus, 
there  is  a  17%  reduction  in  speed  due  to  the  presence  of  the  walls. 


DIMENSIONLESS  FREOUENCY  kdo 

Fig.  11  -  Relative  wave  speed  c/c^  as  a  function  of  dimensionless  frequency 
k^a  for  a  water-filled  tube  of  the  same  cross-sectional  dimensions 
and  material  as  the  G40  calibrator.  Ratio  of  the  cylinder  radii 
is  1.13.  Vertical  arrows  indicate  asymptotes;  c^  is  the 
dilatational  wave  speed,  cy  is  the  bar  speed,  Cq  is  the 
compress ional  wave  speed  in  water. 

The  next  question  concerns  the  calibration  of  a  hot-film  particle  motion 
hydrophone  in  the  calibrator.  In  order  to  increase  the  particle  velocity  for 
a  given  pressure  amplitude,  one  can  reduce  the  propagation  speed  by  choosing  a 
less  rigid  wall  material.  Rubber  might  be  too  weak  to  support  the  water 
column,  but  a  lucite  wall  might  be  realizable.  Let  the  dimensions  be  the  same 
as  for  the  G40.  The  data  for  lucite  are:  density  1200  kg/ra°,  shear  modulus 
0.14x10*°  Pa,  bulk  modulus  0.65x10*°,  dilatational  wave  speed  2650  ra/s.  The 
speed  in  water  for  kja  =0.1  (1035  Hz)  is  reduced  by  a  factor  of  0.30  or  a 


speed  of  440  m/s.  Another  possible  way  to  reduce  the  wave  speed  in  the  tube 
is  to  insert  a  concentric  layer  of  a  suitable  material  inside  the  tube.  The 
parameters  influencing  the  reduction  in  wave  speed  are  the  elastic  constants, 
density,  and  the  ratio  b/a.  The  results  for  a  number  of  materials  are  shown 
in  Table  6.  Only  polyethylene  gives  a  sizable  reduction  in  wave  speed:  a 
factor  of  2.  The  data  for  glass,  lead,  rubber,  and  lucite  were  taken  from 
Ref.  15,  the  data  for  polyethylene  and  polystyrene  from  Ref.  16. 


Table  6  -  Reduction  in  wave  speed  for  various  wall  coatings  in  G40 
calibrator  (Frequency  is  1035  Hz,  a  m  7.0  cm,  b  -  9.8  cm). 


MATERIAL 

G/K 

co/cd 

po/ps 

o 

o 

o 

Glass 

0.641 

0.265 

0.434 

0.99 

Lead 

0. 131 

0.722 

0.088 

0.86 

Lucite 

0.215 

0.559 

0.832 

0.84 

Po lyechy lene 

0.085 

0.760 

0.890 

0.54 

Polystyrene 

0.326 

0.630 

1.065 

0.82 

Rubber  (hard) 

0.200 

0.617 

0.907 

0.79 

Recent  theoretical  and  experimental  work  in  sound  propagation  in  pipes 
with  acoustic  impedance  comparable  to  that  in  water  is  found  in  Refs.  17  and 
18. 


For  comparison  with  the  literature,  Fig.  12  shows  the  dispersion  relation 
in  a  water-filled  tube  in  addition  to  the  curves  for  an  empty  tube,  and  a 
water  column  that  has  a  rigid  boundary  and  one  that  has  a  free  boundary.  The 
material  of  the  wall  is  brass,  as  in  Ref.  19,  the  results  of  which  are 
reproduced  in  Ref.  20.  The  ratio  of  cylinder  radii  was  chosen  accordingly, 
b/a  *  1.13.  The  material  constants  for  brass  are  those  given  in  Ref.  15. 

Ps  »  8500  kg/m3,  G/K  =  0.279,  and  c^  *  4700  m/s.  The  "longitudinal  wave 
speed"  is  given  in  Ref.  19  as  3590  m/s  which  actually  is  closer  to  the  bar 
speed,  3500  m/s.  according  to  Ref  16.  The  results  in  Fig.  12  show  that  the 
wave  speed  approaches  the  bar  speed  (subscript  Y)  for  small  k^a,  in 
disagreement  with  Ref.  19  but  in  accordance  with  Ref.  21. 


Fig.  12  -  Relative  wave  speed  c/c^  as  a  function  of  dlnenslonless 

frequency  kda  for  a  water-filled  brass  tube,  according  to  Ref. 
21,  with  the  wave  speeds  for  the  separate  composing  parts:  tube 
and  water  column.  The  ratio  of  the  cylinder  radii  is  1.13;  cd  is 
the  dllatational  wave  speed,  cY  is  the  bar  speed,  and  cQ  is  the 
compresslonal  speed  in  water.  The  mode  designation  for  the 
fluid  column  follows  the  nomenclature  of  Ref.  6  by  which  the 
first  entry  gives  the  number  of  radial  nodal  lines,  and  the 
second  the  number  of  nodal  circles. 


6.0.  CONCLUSIONS 

The  mathematical  analysis  and  computer  programs  developed  in  this  report 
may  be  used  for  determining  spatial  properties  of  the  acoustic  fields, 
including  wave  speed,  in  Impedance  tube  and  calibrator.  Their  application  to 
the  Interpretation  of  data  from  these  instruments  is  subject  to  the 
restrictive  assumptions  underlying  the  analysis.  First,  in  relating  the 
reflection  coefficient  to  propagation  speed  in  the  impedance  tube  it  is 
assumed  that  the  radial  variation  of  pressure  and  particle  velocity  is  small 
or,  in  other  words,  that  no  other  than  the  zero  mode  of  the  wave  in  the  sample 
needs  to  be  considered.  In  the  second  place,  it  is  assumed  throughout  the 
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computation  of  wave  speeds  that  the  composite  cylinders  are  infinite  in  the 
axial  direction.  This  ignores  the  influence  of  the  unavoidable  termination  in 
a  finite  cylinder  on  the  wave  pattern  in  the  waveguide.  The  experimental 
efforts  to  terminate  the  waveguides  by  rigid  reflectors  or  other  idealized 
devices  can  ipso  facto  be  only  approximately  successful.  It  is  assumed  that 
the  idealization  implied  by  the  model  of  infinite  cylinders  will  be  more 
closely  approached  the  larger  the  ratio  of  length  to  diameter  of  the  sections 
of  the  waveguide. 

Keeping  this  restriction  in  mind,  one  may  conclude  that  it  is  not 
necessarily  true  that  the  wave  speed  measured  by  the  impedance  tube  is  the 
dilatational  speed.  It  is  rather  a  function  of  the  thickness  of  the  fluid- 
filled  gap  between  sample  and  wall,  in  relation  to  the  elastic  properties  of 
the  sample  material  that  determines  the  proper  wave  speed.  As  a  consequence, 
one  may  measure  different  wave  speeds  by  changing  the  radius  of  the  sample,  or 
by  removing  the  fluid  from  the  gap.  Experimental  study  will  have  to  decide  in 
how  far  the  finite  length  of  the  sample  would  interfere  with  the 
interpretation  of  the  data. 
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Coapater  Progri 


In  this  appendix  the  computer  programs  are  described,  and  listings  and 
examples  are  presented.  An  overview  of  the  programs  is  given  in  Table  Al. 


Table  Al  -  Name  and  purpose  of  computer  programs* 


NAME 

TYPE 

INPUT 

OUTPUT 

PURPOSE 

RTUBE 

Complex 

Reflection 

coefficient 

Propagation 

speed 

Impedance 

tube 

WGUIDE 

Options 

1  &  2 

Real 

G/K  and 

co/cd  of 
sample 

Propagation 

speed 

tube 

Speed  in 
sample  in 
impedance 

WGUIDE 

Options  3-7 

Real 

G/K  and 

co/cd  of 
wall 

Propagation 

speed 

Speed  in 
calibrator 

IMPED 

Complex 

G/K  and 

co/cd  of 
sample 

Propagation 

speed 

Speed  in 
sample 

*  G,  shear  modulus;  K,  bulk  modulus;  cQ  speed  in  water;  dilatational  speed 
in  solid. 
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DESCRIPTION  OP  PROGRAM  RTUBE 


The  computer  program  RTUBE  determines  numerical  values  for  the  sound 
speed  and  attenuation  of  a  sample  in  the  acoustic  impedance  tube.  The  input 
to  the  program  consists  of  the  measured  values  of  the  amplitude  and  phase  of 
the  signal  reflected  by  the  sample  and  of  the  signal  reflected  without  a 
sample  in  the  tube,  indicated  by  "standard".  The  program  is  written  in 
FORTRAN  IV  and  should  run  on  any  computer  that  accepts  this  language. 

There  are  two  major  parts  in  RTUBE.  First  a  real-valued  seed  or  starting 
value  for  the  sound  speed  is  computed  from  the  reactive  part  of  the  acoustic 
impedance.*  This  seed  is  entered  into  subroutine  R0819,  which  contains  the 
logic  for  the  complex  root  search.  The  function,  the  roots  of  which  are 
sought,  is  contained  in  subroutine  DS.** 

The  real  root  finder  determines  the  real  part  of  (kd)  corresponding  to 
the  reactive  part  X  of  the  impedance  Z  in  Eq.  (2).  according  to  the  relation 

~T-  *  “CP  /P  )  cot(kd)  .  (Al) 

k  a  so 

o 

Obviously,  there  is  an  infinite  number  of  roots  kd  of  this  equation  and  one 
needs  circumstantial  evidence  about  the  propagation  speed  in  order  to  choose 
the  correct  one.  In  the  iterative  procedure  used  here,  it  is  necessary  to 
evaluate  an  inverse  cotangent  function  which  involves  choosing  the  proper 
branch. 

The  propagation  speed  in  the  water  is  computed  for  a  given  temperature 
and  pressure.  The  density  of  the  fluid  is  1.034  g/cm3,  applicable  to  the 
glycol  water  mixture  used  in  the  experiment.  For  a  different  fill-fluid,  this 
number  should  be  changed  accordingly. 

Important  Variables  in  RTUBE 

Input  Variables 

AMAX  Maximum  frequency  (Hz) 

AMIN  Minimum  frequency  (Hz) 

DENS  Ps>  density  of  material  (g/cra3) 

PRES  p,  pressute  (MPa) 


*  C.M.  Ruggiero,  "Solution  of  Transcendental  and  Algebraic  Equations  with 
Application  to  Wave  Propagation  in  Elastic  Plates,"  NRL  Memorandum  Report 
4449,  November  1981. 

**  P.S.  Dubbelday,  "Complex  Root-Finding  Program  with  Application  to  the 
Dispersion  Relation  of  Waves  Propagating  in  a  Fluid  Plate,"  NRL  Memorandum 
4559,  November  1981. 
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PSA  Phase  of  sample 

PST  Phase  of  standard 

STEP  Step  size 

TEMP  Temperature  (#C) 

THICK  Length  of  sample  (cm) 
VSA  Voltage  of  sample 

VST  Voltage  of  standard 


N  Resistance 

Q  2  kQd ,  kQ  Is  the  wave  number  in  water. 

U  Number  used  in  slope  computation 

V  Absolute  value  of  reflection  coefficient 

W  Wave  speed  in  water  corrected  for  ambient  temperature 

and  pressure 

Y  cot(kd) 
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Output  variables 


AMGA  Attenuation  (Np/m) 

FREQ  Frequency  (Hz) 

SS  Sound  speed  (m/s) 

Major  Computation  Blocks  in  RTUBE ,  10819 

Descriptions  of  the  major  computation  blocks  in  RTUBE  are  as  follows 

Line  Number 


Computation  Block 

From 

To 

Input  section 

2 

25 

Slope  computation 

41 

42 

Iterative  rootfinder 

43 

59 

Output  section  (real  roots) 

64 

66 

Descriptions  of  the  major  computation  blocks  in  R0819  are  as  follows 


Line 

Number 

Computation  Block 

From 

To 

Set  vertical  pair 

31 

32 

Test  for  initial  location  of 

39 

40 

vertical  pair 

Determination  and  movement  of 

41 

49 

vertical  test  pairs 

52 

53 

Check  for  iterations 

62 

63 

Check  for  sign  change  in 

67 

vertical  pair 

Set  horizontal  pair 

68 

Determination  and  movement  of 

74 

75 

horizontal  test  pairs 

86 

88 

Check  for  sign  change  in 

102 

horizontal  pair 

Check  for  desired  accuracy 

III 

112 

Reduce  step  size 

114 

34 


Horizontal  vertical  error 
(real  &  imaginary) 

120 

121 

Advancing  y  value 

124 

Output  section 

127 

131 

I 
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MM Mft 

*"  *  *  -  *  *  *  *  ■  ’  '  . 

m 

c 

c 

RTUBE 

c 

THIS  PROGRAM  COMPUTES  THE  SOUND  SPEED  AND  ATTENUATION 

c 

OF  THE  IMPEDANCE  TUBE  WHEN  THE  VOLTAGE  OF  THE  SAMPLE  AND 

c 

THE  VOLTAGE  OF  THE  STANDARD  >  THE  PHASE  OF  THE  SAMPLE  AND  THE 

c 

PHASE  OF  THE  STANDARD  ARE  MEASURED. 

c 

SUBROUTINES  USED! 

■1 

c 

c 

R0819 

DS 

■ 

0001 

REAL  M»N»I.J.L»DAT!5>9) 

0002 

WRITE ! 5 » 10 > 

Br/  CS 

000  3 

10 

FORMAT ( '  ENTER  PRESSURE ! MPA >. LENGTH ! CM >» TEMP ! C >» DENSITY ! G )  ') 

0004 

READ(5»20)PRES.THICK»TEMPiDENS 

a 

0005 

20 

F0RMAT<4F10.0) 

0006 

THICK=THICK/100 

0007 

ZNOT  =  2.5 

r% 

0009 

PI  =  0 

0009 

ANS 1 =- 1 0000 

0010 

K  =  1 

1  .  • 

0011 

WRITE!5»21) 

0012 

21 

FORMAT!/.'  ENTER  FREQUENCY  < MIN . MAX . STEP  )  ') 

0013 

READ (5. 23) AMIN rAMAX. STEP 

** 

0014 

23 

FORMAT !  3  F 1 0 . 0  ) 

k: 

0015 

0PEN!UNIT=1 , NAME3 'DATXY.DAT' » T YPE= ' NEW ' » FORM3 ' UNFORMATTED ' ) 

0016 

J1  =  0 

0017 

DO  500  FPEQ=AMIN, AMAX.STEP 

- 

0019 

PI  =  0 

' 

0019 

J1=J1+1 

0020 

WRITE!5,25)FREQ 

0021 

25 

FORMAT!//.'  FREQUENCY3  '.FIO.l./) 

0022 

29 

WR I TE ( 5  »  30  ) 

0023 

30 

FORMAT!'  ENTER  VOLT.  STAND.. VOLT.  SAMP.  PHASE  STAND.. PHASE  SAMP '  ) 

0024 

READ!5.40,ERR=29)  VST . USA , PST . PSA 

0025 

40 

FORMAT  <  4F10 . 0 ) 

0026 

V=VSA/VST 

0027 

W  = 1 460. 0+ 1 0. 0*PRES/6. 96+ 2. 8* ! TEMP-24. 0)-. 037*! TEMP-24.0) **2.0 

oo2e 

0=4*3. 14159*THICK*FREG/W 

0029 

E=!PSA-PST)/57.2958 

m 

0030 

G  =  E-G 

au 

0031 

IF  (ABS!G)-2*3. 14159  .LE.  0)  GO  TO  50 

0  0  3  2 

IF  (G  .LE.  0 ) GO  TO  60 

0033 

G=G-2*3. 14159 

0  0  3  4 

GOTO  50 

0035 

60 

G=G+2*3. 14159 

0036 

50 

N= ( 1 - V**2 ) / ( 1+V**2-2*V*C0S ! b ) ) 

0  03  7 

M  =  2*V*SIN(G>/! 1+V**2-2*V*C0S!G)  ) 

• 

>0  38 

1.0  SHOULD  BE  CHANGED  TO  1.034  FOR  GLYCOL-WATER 

A=DENS*2*3. 14159*FREG*THICK/(W*1 .0) 

?  v>  3  g 

J= . 0001 

00  40 

U=. 00001 

00  4  1 

90 

I=-M*!ZNOT-U)+M*(ZNOT+U) 

0042 

X=A*!COS(ZNQT-U)/SIN(ZNOT-U> ) - A* ( COS ( ZNOT +U ) /S I N ! ZNOT +U ) ) 

>  V  **  T 

ASSIGN  95  TO  JOB 

• 

• 

'• 

36 

0044 

IF  ( ABS  < I )  .LE.  ABS ( X  > 1  ASSIGN  100  TO  JOB 

0045 

GO  TO  JOB 

0046 

95 

Y  =  A*C0S(ZN0T  +  F'I*3 .14159  ) /SIN ( ZN0T+PI*3 . 14 1 ! 

0047 

Z=Y/<-M> 

0048 

I F ( K  .GT.  1000 (ASSIGN  100  TO  JOB 

004? 

GO  TO  110 

0050 

100 

H=-M*(ZN0T+FI*3. 14159  >/A 

0051 

Z=ACOS(H/SORT< 1+H**2>  > 

0052 

I F  (  K  .GT.  1 000  >  ASSIGN  95  TO  JOB 

0053 

110 

DIFF=ZNOT-Z 

0054 

IF ( ABS <  D IFF )  .LE.  J)  GO  TO  200 

0055 

ZNOT=Z 

0056 

K  =  K+  1 

0057 

IF  <  K  .GT.  1001 >ZN0T=ANS1 

0058 

IF  (K  .GT.  1001 )K=1 

0059 

GOTO  JOB 

0060 

200 

ANS=Z+PI*3. 14159 

0061 

IF  (ANSI  .GT.  CANS  >>PI=1 

0062 

IF ( ANSI  .GT.  ANS>ZN0T=ANS1 

0063 

IF (ANSI  .GT.  ( ANS  >  >  GOTO  JOB 

0064 

URITE(5.210)ANS.M.N 

0065 

210 

FORMAT (  '  REAL  SEEDi  REACT ANCE  = i  RESISTANCE* 

0066 

URITE< 1 )FREQiH»N. A. ANS 

0067 

ANS1=ANS 

0068 

500 

CONTINUE 

0069 

CLOSE <  UN  I T* 1 ) 

0070 

CALL  R0819<J1 (THICK) 

0071 

END 

r-'; 


rrrr'rrm 

Wi'V 


K. 


1 


ta 


0001 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE  R0819 < J 1 » THICK > 

THE  PURPOSE  OF  THIS  SUBROUTINE  IS  TO  DETERMINE  THE 
COMPLEX  ROOTS  OF  FUNCTION  DS 

DESCRIPTION  OF  PARAMETERS 

J1  -THE  NUMBER  OF  FREQUENCIES  TO  BE  COMPUTED 
THICK  -THE  THICKNESS  OF  THE  ELEMENTS 


FUNCTION  REQUIRED 

Z  -  APPROXIMATE  ROOT 


USAGE 


note: 


DS(Z) 


ALL  DATA  IS  PASSED  IN  A  COMMON  STATEMENT^ 


0002 

COMPLEX  All » A22»  SI  1 » S22 > DELA » DELS » Z » F  » DSP  > AI 

lZlfZ2.Z3.Z4,Fl»F2»F3.F4»CSTH»FF>DSF2.AKAP»AI 

0003 

COMPLEX  DS»CSQRTfCMPLX.CCOS.CSIN» CHECK  .TEST 

0004 

REAL  LOOK 

0005 

COMMON  R'HOrAKSD.A 

0006 

OPEN ( UNI T= 1 .NAME='DATXY .DAT ' , TYPE= ' OLD ' t FORM 

0007 

K  =  1 

0008 

ANR  =  8 

0009 

TOL=  .  0 1 

0010 

REDF=10 

0011 

DELZ= • 0 1 

0012 

KAM=10000 

0013 

DO  10  I=1»J1 

0014 

READ  < 1 >  FREQ  r AKSD  »RHM » A  >  ZNOT 

0015 

ZNOT=ZNOT/THICK 

0016 

RH0  =  0 

0017 

AI -CMPLX  <  0  *  »  1  .  ) 

0018 

NR  =  I  NT ( ANR ) 

0019 

NRV  =  0 

c 

COMPUTE  STEPSIZE  DELZ 

0020 

ZUAR=ZNOT 

0021 

Z=ZVAR+DELZ 

0022 

DFDZ=OS<Z)/DELZ 

0023 

RHO=RHO+RHM/ANR 

002  4 

NRU=NRU+1 

0025 

DZ=-DS ( ZVAR ) /DFDZ 

0026 

X  =  ABS ( REAL  <  DZ )  ) 

0027 

Y  =  ABS< AIMAG(DZ)  ) 

0023 

DEL=AMAX1<X»Y>*1 

00  29 

DELP=DEL 

0030 

70 

DEL=DELP 

0031 

64 

Z1=ZVAR-AI*DEL 

0032 

Z2=ZUAR+AI*DEL 

0033 

65 

F1=DS(Z1) 

0  0  3  4 

F2=DS(Z2) 

0  0  3  5 

NR  =  0 

00  3  o 

KL  =  0 

38 


.1 


0037 

KUP=0 

0038 

KDOWN-O 

0039 

SIGN-AIMAG<F1)*REAL(F2)-REAL<F1)*AIMAG<F2) 

0040 

IF(SIGN)7l. 72.73 

0041 

71 

FH1--1. 

0042 

GO  TO  74 

0043 

72 

WRITE(S, 302)  SIGN 

0044 

302 

FORMAT  ('  LEFT-RIGHT  SIGN  *  '.E12.5) 

0043 

GO  TO  92 

0046 

73 

FH1-1. 

0047 

74 

DELH-FHltDEL 

0048 

75 

Z1*Z1+PELH 

0049 

Z2-Z2+DELH 

0050 

F1«D8(Z1> 

0051 

F2*D8(Z2) 

0052 

SIGN-AINAG(F1)*REAL(F2>-AIMAG(F2)*REAL(F1> 

0033 

IF  (SIGN)  81 1 82 1 83 

0054 

81 

FH2--1. 

0055 

KL«KL+1 

0056 

IF  (KAM-KL)  99.99.84 

0057 

99 

URITE<5.310) 

0058 

310 

FORMAT  ('  EXIT  LEFT') 

0059 

GO  TO  100 

0060 

82 

GO  TO  92 

0061 

83 

FH2  -  1. 

0062 

KR*KR+1 

0063 

IF (KAM-KR)  98.98,84 

0064 

98 

WRITE  (3.311) 

0065 

311 

FORMAT <  '  EXIT  RIGHT' > 

0066 

GO  TO  100 

0067 

84 

IF  (FH1*FH2>  92,75,75 

0068 

92 

Z4=(Z14Z2>/2.+DEL 

D 

WRITE  (5,886)  KR.KL 

DB86 

FORMAT  ('  KR» ' » 14 , '  KL-'.I4) 

0069 

KR-0 

0070 

KL*0 

0071 

Z3=(Zl+Z2>/2. -DEL 

0072 

F3»DS(Z3) 

0073 

F4*DS(Z4) 

0074 

SIGN=AIMAG (F3)*REAL(F4) -AIMAG (F4)*REAL(F3) 

0073 

IF  (SIGN)  101,102,103 

0076 

101 

F\>1  =  1. 

0077 

GO  TO  104 

0078 

102 

WRITE  (3,304)  SIGN 

0079 

GO  TO  122 

0080 

103 

FV1  =  -1  . 

0081 

104 

DELV«FV1*DEL 

0082 

105 

Z4»Z4+AI*DELV 

0083 

Z3=Z3+AI*DELV 

0084 

F3=DS(Z3) 

0085 

F4*DS(Z4) 

0086 

SIGN-AIMAG(F3)*REAL(F4)-AINAG<F4)*REAL(F3> 

0087 

304 

FORMAT ( '  UP-DOWN  SIGN  =',E12.3> 

0088 

IF  (SIGN)  111,112,113 

0089 

111 

FV2  -1. 

KUP-KUP+1 

IF<  KAH-KUP )  97.97.114 
WRITE  (5.305 > 

FORMAT  ( '  EXIT  UP'  ) 

GO  TO  100 
GO  TO  122 
FV2  —  1 . 

KDOUN-KDOWN+1 
IF  (KAM-KDOWN)  96.96.114 
WRITE ( S  >  306 ) 

FORMAT  < '  EXIT  DOWN') 

GO  TO  100 

IF  <F01*FV2)  122.105.105 
REGA=REAL(Z1 > 

WRITE  <  5 . 888 )  KUP.KDOWN 

FORMAT  <'  KUP“ '  «  14 «  '  KD0MN='.I4> 

AMGA'AIMAG ( Z3 ) 

KUP=>0 

KD0UN=0 

IF  ( AMGA >  133.134.133 
WRITE  <5.233) 

FORMAT  <  '  AMGA  IS  ZERO') 

GO  TO  131 
AMGA-ABS(AMGA) 

IF (DEL/AMGA-TOL >132.132. 131 

CONTINUE 

DEL=DEL/REDF 

Zl=<Z3+Z4)/2.-AI*DEL 

Z2=<Z3+Z4)/2.+AI*DEL 

GO  TO  65 

REGA=REAL  <Z1>  -DELH/2. 

ANGARA I MAG  (Z3)-DELU/2. 

DELH=ABS < DELH/2  .  ) 

DELU=*ABS<DELU/2. ) 

IF (NR-NRU >143.145. 146 

ZVAR-REQA4AMGABAI 

RHO-RHO+RHM/ANR 

NRY=NRU+1 

GO  TO  70 

WRITE  (5.173 ) FREQ 

FORMAT!/.'  SOUND  SPEED  AND  ATTENUATION  FOR  FREQ 
SS=<2*3. 14159*FREQ)/REGA 
WRITE  <5»174)SS»ABS<AMGA) 

F0RMAT(2E15.5) 

CONTINUE 
CLOSE ( UN I T  =  1  ) 


0001 

FUNCTION  DS(Z> 

c 

c 

THIS  IS  A  SPECIFIC  FUNCTION  TO  BE  CALCULATED  UITH 

c 

p 

COMPLEX  ROOTFINDER.  THE  COMPLEX  ARGUMENT  IS  Z. 

L 

c 

NOTE 

c 

ALL  DATA  IS  PASSED  IN  COMMON  STATEMENT 

0002 

COMPLEX  DSiZ»All.A12pA22»AC0TZ.C0TZ»SINZ»C0SZ»CC0S»CSINrSH.CMPLX 

0003 

COMMON  RHOf AKSDiTHICK.A 

0004 

13 

F0RMAT<4F10.3»2E15.5> 

0005 

AI=CMPLX<0. 0»1.0) 

0006 

A12=  -THICK*Z 

0007 

All=CMPLX(-AKSDiRHO) 

0008 

A22*A12*A11 

0009 

SI 1 =THICK*Z 

0010 

C0SZ=CC0S(S11) 

0011 

S I NZ'CSI N ( S 1 1  ) 

0012 

COTZ=COSZ/S INZ 

0013 

ACOTZ=A*COTZ 

0014 

DS*AC0TZ+A22 

0015 

RETURN 

0016 

END 

EXAMPLE 


RUN  RTUBE  <CR> 

ENTER  PRESSURE<MPA)»LENGTH<CM),TEMP<C>»DENSITY<G) 
.61,12.09,23.9,1.13  <CR> 

ENTER  FREQUENCY  < MIN . MAX » STEP > 

4000* 11000. 1000  <CR> 


FREQUENCY®  4000.0 

ENTER  VOLT  .  STAND. .VOLT.  SAMP.  PHASE  ST AND .. PHASE  SAMP 
.084. .077.269.91,387.41  <CR> 

REAL  SEED.  REACTANCE®,  RESISTANCE®  2.22014  0.77738 


FREQUENCY®  SOOO .  0 

ENTER  VOLT.  STAND. .VOLT.  SAMP,  PHASE  STAND. .PHASE  SAMP 
.127,  .112, 216. 6.237, 9_A  <CR> 

REAL  SEED,  REACTANCE®,  RESISTANCE®  2.76703  2.61237 


FREQUENCY®  6000.0 

ENTER  VOLT.  STAND., VOLT.  SAMP,  PHASE  STAND., PHASE  SAMP 
■183. .148,176.40.192.52  <CR> 

REAL  SEED.  REACTANCE®.  RESISTANCE®  3.33705  -4.64150 


FREQUENCY®  7000.0 

ENTER  VOLT.  STAND., VOLT.  SAMP,  PHASE  STAND., PHASE  SAMP 
.231 . .173,136.86,145.43  <CR> 

REAL  SEED.  REACTANCE®,  RESISTANCE®  3.77241  -1.44402 


FREQUENCY®  8000.0 

ENTER  VOLT.  STAND., VOLT.  SAMP,  PHASE  STAND., PHASE  SAMP 
■  263 . .  177 ,101.86,11  1.03  <CR> 

REAL  SEED,  REACTANCE®,  RESISTANCE®  4.28325  -0.48575 


FREQUENCY®  9000.0 

ENTER  VOLT.  STAND ., VOLT .  SAMP,  PHASE  STAND. .PHASE  SAMP 
. 308, . 190,71 .89,85.73  <CR> 

REAL  SEED,  REACTANCE®,  RESISTANCE®  4.79123  0.08435 


FREQUENCY®  10000.0 

ENTER  VOLT.  STAND., VOLT.  SAMP,  PHASE  STAND., PHASE  SAMP 
.338, .220,40. 17,53.37  <CR> 

REAL  SEED,  REACTANCE®.  RESISTANCE®  5.24966  0.64499 


FREQUENCY®  11000.0 

ENTER  VOLT.  STAND., VOLT.  SAMP,  PHASE  STAND., PHASE  SAMP 
.339.  ,242’35?-.34.°.79  <CR> 

REAL  SEED,  REACTANCE®.  RESISTANCE®  5.72622  1.75329 


0.06983 


0.50517 


4.21691 


0.46431 


0.24388 


0.23873 


0.30552 


0.75410 


4000 


SOUND  3  F'  E  E  D  AN  If  A  I  TENUATZCN  FOP.  FREQ  =! 
0.13677E+04  0.45/41E+00 

SOUND  SPEED  AND  ATTENUATION  FOP  FREO  =:  5000. 

0. 13656E+04  0.60050E+00 

SOUND  SPEED  AND  ATTENUATION  FOR  FREQ  =:  6000. 

0.13943E+04  o.asoosE^oo 

SOUND  SPEED  AND  ATTENUATION  FOR  FREQ  7000. 

0.14196EtD4  C. ;v334F+01 

SOUND  SPEED  AM  ATTFNJh': ION  FOR  FREO  =1  8000. 

o .  1405  5  etc  a 

SOUND  SPEED  AND  aT’En'JATIQN  FOR  FREQ  = :  9000. 

0  .  1  4290E-'- 1  '  J.1VC29EP01 

SOUND  SPEED  AND  AT'ENUATION  FOR  FREQ  =1  10000. 

0.14401E  +  04  <:■ .  1  9  o  l  9  E  +  0 1 

SOUND  SPEED  AND  AT’ENUATION  FOP  FREQ  = ’.  11000. 

0.14439ER04  0 . i57?9Et01 


1  >  All  underlined  r  irtlons  are  user  supplied. 

2)  'CH>  indicates  ’RETURN' 


DESCRIPTION  OP  PROGRAM  WGUIDE 


The  program  WGUIDE  computes  the  speed  of  axially  symmetric  waves 
propagating  along  the  axial  direction  in  composite  Infinite  cylindrical  wave¬ 
guides,  given  the  values  of  the  ratio  of  shear  modulus  and  bulk  modulus  of  the 
solid  and  the  ratio  of  the  propagation  speed  in  the  fluid  to  the  dilatatlonal 

wave  speed  in  the  solid.  To  facilitate  tracing  the  complete  graph  for  the 

dispersion  relation,  Including  higher  modes,  two  methods  of  scanning  are 
possible.  In  the  first  one,  the  dimensionless  frequency  k^a  (abscissa)  is 
fixed  and  the  program  searches  for  all  relative  wave  speeds  (ordinate)  between 
two  limits.  In  the  other  method,  the  program  steps  through  increasing  or 
decreasing  values  for  the  dimensionless  frequency  and  finds  the  smallest  value 
for  the  relative  wave  speed  that  falls  between  limits,  chosen  such  that  the 
desired  branch  is  bracketed.  It  is  also  possible  to  determine  the  vertical 
asymptotes.  The  application  of  the  program  is  twofold:  firstly,  the  program 
computes  the  wave  speed  in  the  sample  Inside  an  impedance  tube  which  is 

assumed  to  be  rigid  and  separated  from  the  tube  wall  by  a  concentric  layer  of 

fluid;  secondly,  one  may  apply  the  program  to  assess  the  influence  of 
nonrigidity  of  the  wall  for  inertial  calibrators. 

The  main  program  calls  subroutine  FUNC  for  Option  numbers  1  through  5. 
FUNC,  in  turn,  calls  subroutine  DI3  which  contains  the  impedance  tube 
dispersion  elation,  for  Options  1  and  2.  Number  1  assumes  that  the  impedance 
tube  wall  is  rigid,  number  2  (rather  unrealistically)  assumes  a  free  outer 
boundary.  For  Option  numbers  3,  4,  and  5,  FUNC  calls  subroutine  CALIB,  which 
contains  the  dispersion  relations  for  the  inertial  calibrator  arrangement. 
Number  3  is  for  an  elastic  tube,  fluid-filled  or  empty;  number  4  is  for  a 
rigid  tube  with  inner  coating  cemented  to  the  wall;  and  number  5  is  for  a 
rigid  tube  with  inner  coating  free  to  move  tangentially  with  respect  to  the 
wall.  If  the  Option  number  is  6  or  7,  the  main  program  calls  subroutine 
FLUID;  for  number  6,  it  computes  the  wave  speed  for  a  fluid  tube  with  rigid 
boundary.  Note  that  the  zero  order  is  the  compresslonal  wave  speed  in  the 
fluid.  For  the  number  7,  the  wave  speed  in  a  fluid  with  free  boundary  (in 
this  case  there  is  a  low-frequency  cut-off)  is  calculated. 


The  program  computes  the  location  of  vertical  asymptotes  by  entering  a 
negative  value  of  AKDA.  In  Option  1,  entering  a  negative  value  fo  RHM  will 
prompt  the  computation  of  the  wave  speed  in  a  solid  cylinder  in  vacuum 
(indicated  by  "VACUUM  TUBE"  in  the  listing  and  printout).  In  Option  3, 
setting  RHM  equal  to  0  will  result  in  the  computation  of  the  wave  speed  in  an 
empty  cylindrical  shell,  indicated  by  "EMPTY  TUBE"  in  print-out. 


The  Bessel  function  subroutines  are  documented  in:  System/360  Scientific 
subroutine  Package  Version  II  (c).  International  Business  Machines  Corporation 
1966,  1967,  1968,  White  Plains,  NY,  Section  -  Special  Operations  and 
Functions,  363  -  367  (copies  of  these  subroutines  will  not  be  found 
in  this  report;  but  the  user  may  implement  the  needed  Bessel  function  routines 
from  the  program  library).  Subroutine  DET  computes  the  determinant  value  of  a 
matrix.  The  various  results  are  stored  by  calling  subroutine  FILE,  which  is 
prompted  by  choosing  Option  number  9.  The  values  of  0^  “  c^/c., 

SR  ■  c^/cj,  0Q  ■  Cq/Cjj,  and  the  asymptotes  are  also  stored  in  a 


0 


cjc  , , 
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file  separate  from  the  dispersion  results  in  order  to  discriminate  between 
them.  They  ran,  however,  appear  on  the  same  graph  by  running  a  plotting 
program  named  PLOTTER.* 


FORTRAN  Variable  Names  for  WED IDE 


Input  Variables 


KB:  Option  number  determining  mode  of  operation 

*  1  Sample  in  fluid  tube  with  rigid  boundary 
=  2  Sample  in  fluid  tube  with  free  boundary 

-  3  Fluid-filled  elastic  tube  or  empty  tube 

=  4  Fluid-filled  coated  rigid  tube  (cemented) 

=  5  Fluid-filled  coated  rigid  tube  (not  cemented) 

*  6  Fluid  tube  with  rigid  boundary 
=7  Fluid  tube  with  free  boundary 

GOK  (G/K)  ratio  of  shear  modulus  to  bulk  modulus 

COCD  (Cq/c^)  ratio  wave  speed  in  fluid  to  dilatational 

wave  speed  in  solid 

REM  (p  / p  )  ratio  of  fluid  density  to  density  of  solid;  negative 

value  for  solid  cylinder  in  vacuum;  zero  for  empty  calibrator 

AKDA  (k^a)  dimensionless  wave  number  of  dilatational  waves 

,  u>a 

k.a  =  — 

d  °d 

Bi  or  BAM  Minimum  (k^a  if  asymptote  option,  3  otherwise) 

B2  Maximum  (k^i  if  asymptote  option,  3  otherwise) 

DB  or  DEL  Step  size  of  (k^a  if  asymptote  option,  3  otherwise) 

KO:  Return  option  number 

=  i  program  returns  to  entering  input  option  number 

-  2  program  returns  to  entering  GOK 

=  3  program  returns  to  entering  COCD 

=  4  program  returns  to  entering  RHM 

-  5  program  returns  to  entering  MIN,  MAX,  STEP  (Y  value) 

-■  o  program  returns  to  entering  BDA 

=  7  urogram  returns  to  entering  AKDA 

=  o  program  returns  to  entering  MIN,  MAX,  STEP  (X  value) 

=  9  stores  dat  i 
10  program  exits 


*  C.M.  Ruggiero  and  R.W.  Anderson,  "General  Purpose  High-Resolution  Plotting 
Package  for  Tektronix  4662  Plotter  and  Compatible  CRT  Terminals,"  NRL 
Memorandum  Report  4333,  August  1131. 


A 

B 

Dl  and 
TOL 

Output 

AKF 

ERR 

AKDA 


Array  that  stores  k^a 
Array  that  stores  B 
D2,  Values  of  the  determinant 
Tolerance 

Variables 
Asymptote  at  k^a 
Error  (DEL/2) 
kda 


BAF 


B  ■  c/c. ,  final  value 


Major  Coaputation  Blocks  in  WGOIDR 

Descriptions  of  the  major  computation  blocks  in  WGUIDE  are  as  follows 

Line  Number 

Computation  Block  From  To 

4  55 

totes  61  78 
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SOURCE  LANGUAGE  LISTING 


0001 

0002 

0003 

0004 

0005 


0006 

0007 

0008 

0009 

0010 

0011 

0012 

0013 

0014 

0015 

0016 

0017 

0018 

0019 

0020 

0021 

0022 

0023 

0024 

0025 

0026 

0027 

0028 

0029 


0030 

0031 

0032 

0033 

0034 

0035 


C  WGUIDE 

C  WRITTEN  BY  P.  DUBBELDAY  AND  T.  RUGGIERO 

C 

C  THE  PURPOSE  OP  THIS  PROGRAM  IS  TO  CALCULATE  THE  SPEED 

C  OF  AXIAL  WAVES  IN  CYLINDRICAL  WAVE  GUIDES. 

C 

C  C  SUBROUTINES:  BESJ.BESY.BESK.DET.FILE.FUNC.INUE 

C  I0.DI3. CALIB.FLUID 

C 

REAL  A (500  >  » B  ( 500 ) 

COMMON  AKA »AKDA» BAMS. ACC » COCD » BDA . RHM » KB . AA 
K9  =  0 

110  WRITE  <  5 . 120 ) 

120  FORMAT ( '  ENTER  #  CORRESPONDING  TO  OPTION  DESIRED'. 

1/.'<1)  SAMPLE  IN  FLUID  TUBE  WITH  RIGID  BOUNDARY'. 

2/. '(2)  SAMPLE  IN  FLUID  TUBE  WITH  FREE  BOUNPARY'. 

3/. '(3)  FLUID  FILLED  ELASTIC  TUBE.  OR  EMPTY  TUBE'. 

4/. '(4)  FLUID  FILLED  COATED  RIGID  TUBE  (CEMENTED)'. 

5/. '(5)  FLUID  FILLED  COATED  RIGID  TUBE  (NOT  CEMENTED)'. 

6/. '(6)  FLUID  TUBE  WITH  RIGID  BOUNDARY'. 

7/. '(7)  FLUID  TUBE  WITH  FREE  BOUNDARY') 

READ  (5.640)  KB 

IF ( KB  .GE.  6 ) CALL  FLUID ( KB ) 

TOL* . 00 1 
ACC* . 0 1 

130  WRITE  (5.140) 

140  FORMAT (/. '$ENTER  RATIO  OF  SHEAR  TO  BULK  MODULUS  (G/K)  ') 

READ  (5.170)  GOK 
I F ( K9  . EO .  1 ) GO  TO  300 
150  WRITE  (5.160) 

160  FORMAT  ( ' *ENTER  RATIO  OF  FLUID  TO. SOLID  SPEED  (CO/CD)  ') 

READ (5.170)  COCD 
170  FORMAT  (F15.0) 

IF ( K9  .EO.  1 ) GO  TO  300 
130  UR ITE ( 5  » 1 90 ) 

190  FORMAT ( / »  '  ENTER  RATIO  OF  FLUID  TO  SOLID  DENSITY  '. 

1/.'$  (FOR  EMPTY  TUBE  SET  RHM*O.VACUUM  TUBE  SET  RHM*NEG  VALUE) 

READ(5.170)RHM 

I F ( K9  .Ed.  1 > GO  TO  300 

BAMS*SQRT(G0K/(1 .+4.*G0K/3. ) ) 

K0  =  7 

BY=BAMS*3./SQRT(3.+G0K) 

P0I=(3.-2.*G0K)/(6.+2.*G0K> 

BRE  = (0.87+1. 12*P0I >*BAMS/( 1 . +POI ) 

WRITE(5»200)BAMS»BY.BRE 

200  FORMAT ( / .  '  SHEAR=  '.EB.3.'  YOUNG*  '.E8.3.'  RAYLEIGH*  '.E8.3) 

C 

C  IN  THE  FIRST  3  STORAGE  LOCATIONS  ARE  POINTS  OF  THE  YOUNGS 

C  MODULUS.  IN  THE  4TH  LOCATION  CO/CD  IS  STORED  AT  .1 

A( 1 >*0.5 
B ( 1 ) *COCD 
A  (  2  >  =  1  . 

B ( 2 ) =COCD 
A  <  3 ) =2 . 

B ( 3 ) *COCD 
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0036  A  (  4  )  =  •  1 

0037  8 ( 4  )  =BY 


C 

C 

INPUT  SECTION 

0038 

L 

N  =  4 

003? 

O 

Cl 

UR  I TE ( 5  *  290 ) 

0040 

RE AD  <  5  >  260 ) ANDA 

0041 

IFIK9.EQ.  11G0T0  300 

0042 

IF  <  A  K  D  A  .GE.  0  >  GO  TO  230 

0043 

UR  I TE ( 5  «  220 ) 

0044 

220 

FORMAT (/ . ' $ENTER  MIN.MAX.STEP  (HORIZONTAL  SCAN)  ') 

0045 

GO  TO  250 

0046 

230 

WRITE  (5.240) 

0047 

240 

FORMAT  ( / * ' CENTER  MIN.MAX.STEP  (VERTICAL  SCAN)  ') 

0048 

250 

READ  (5.260)  B1.B2.DB 

0049 

260 

FORMAT (3E15.0) 

0050 

I F ( K9  .EQ.  i  ) GOTO  300 

0051 

270 

WRITE  (5.280) 

0052 

280 

FORMAT  ( / . ' *ENTER  RATIO  OF  TUBE  TO  SAMPLE  RADIUS  (B/A)  ') 

0053 

READ  (5.170)  SB A 

0054 

TYPE  *.  '  ' 

0055 

290 

FORMAT (/>' *ENTER  X  VALUE  (SET  KDA  NEGATIVE  FOR  ASYMPTOTES 

0056 

300 

IF  (ANDA)  310.440.440 

0057 

310 

K0  =  8 

0058 

A  N  D  A  =  B  1 

0059 

320 

DEL  =  DB 

0060 

N  =  N+  1 

C 

c 

FUNC  IS  A  SUBROUTINE  THAT  CALLS  THE  CORRECT  SUBROUTINE  TO 

c 

EVALUATE  THE  DETERMINANT.  KB  IS  THE  OPTION  DESIRED. 

0061 

c 

CALL  FUNC  <  —  1.KB.D1) 

0062 

330 

ANDA*AKDA  +DEL 

0063 

CALL  DECLAR 

0064 

340 

CALL  FUNC ( - 1 . KB . D2 ) 

0065 

350 

IF<D1*D2)  370.540.390 

0066 

360 

GO  TO  620 

006  7 

370 

IF  (ABS(DEL/AKDA)-TOL)  420.380.380 

.’0  6  9 

330 

akda=akda-del 

0069 

DEL  =  DEL  7 10 . 

0  0  7  0 

GO  TO  330 

C 

CHECKS  IF  <  VALUE  EXCEEDS  MAXIMUM  CHOSEN 

00  7  1 

390 

IF  (AKDA-B2/  330.400.400 

0  0  7  2 

400 

W r  !  T  E  ’  5 . 4  1  0  .■ 

TO7  3 

4  1  0 

FORMAT  NO  MORE  ASYMPTOTES  ') 

J  0  7  4 

SC  "0  620 

0  0  75 

4  20 

E  f  -  R  ”  2  £  i-  '  2  . 

0  0  7  6 

•VF  -- AKDA-ERR 

0  0  7  7 

write  <5*430)  A  K  F  .ERR 

0078 

4  1  0 

FCF^AI  .  ASYMPTOTE  AT  X  =  '.F10.5.'  +/-  '.E?.3> 

•  i'[  AL  Y  -*F'T0TE  ARE  FOUNt  AND  POINTS  STORED  AT  Y  =  l. 25>1.5 


0079 

A  ( N )  =  AKF 

00B0 

B(N)=1 .25 

0081 

N=N+ 1 

0082 

A  ( N ) *AKF 

0083 

B ( N )  =  1 . 5 

0084 

GO  TO  320 

0085 

440 

BAM*=B1 

0086 

450 

DEL=DB 

0087 

N  =  N+1 

0088 

CALL  FUNC<  BAM  »  KB  »  D1 ) 

0089 

460 

8AM=BAM+DEL 

0090 

CALL  DECLAR 

0091 

470 

CALL  FUNC(BAM,KB,D2> 

0092 

480 

IF  (D1*D2>  490.540.510 

0093 

490 

IF  ( ABS<DEL/BAM)-TOL)  560.500.500 

0094 

500 

8AM*BAM-DEL 

0095 

DEL=DEL/10 . 

0096 

GO  TO  460 

0097 

510 

IF  (BAM-B2)  460.520.520 

0098 

520 

WRITE  <5. 530) 

0099 

530 

FORMAT  ('  NO  MORE  ROOTS') 

0100 

GO  TO  620 

C 

c 

PRINTED  WHEN  EITHER  DETERMINANT  IS 

EQUALTO 

ZERO 

0101 

540 

WRITE  (5.550) 

0102 

550 

FORMAT  ('  D14D2-ZER0') 

0103 

GO  TO  620 

0104 

560 

ERR=DEL/2. 

0105 

BAF=BAM-ERR 

0106 

570 

IF ( KO  .NE.  7 )G0  TO  600 

0107 

580 

WRITE  (5,590)  BAF.ERR 

0108 

590 

FORMAT ( '  Y  *  '  » F 10 . 5 . '  ERROR  * 

E9.2  ) 

0109 

A ( N ) =AKD A 

0110 

B ( N ) =BAF 

0111 

GO  TO  450 

0112 

600 

WRITE  (5,610)  AKDA, BAF.ERR 

0113 

610 

FORMAT ( '  X  =  '  .  F10 . 5 » '  Y  »  ' , F10 • 5 , 

'  ERROR 

*  '  1 1 

0114 

A ( N ) =AKDA 

0115 

B(N)*BAF 

0116 

AKDA=AKDA+DKHA 

0117 

IF  ( (AKM-AKDA)*DKDA)  620.440.440 

0118 

620 

K9  =  l 

C 

TO  RECYCLE  BACK  OR  EXIT  CHOSE  OPTION  * 

0119 

WR ITE ( 5  »  630 ) 

0120 

630 

FORMAT ( / » '  ENTER  *  CORRESPONDING  TO 

CHANGE 

/ 

1/,'  1)  OPTION' 

2/,'  2)  SHEAR  MODULUS/BULK  MODULUS 

(G/K) 

3/,'  3)  SPEED  IN  FLUID/SPEED  IN  SOLID  (CO/CD! ' 

4/,'  4)  FLUID  DENSITY/SOLID  DENSITY 
5/,'  5)  MIN  MAX.  STEP  <Y  VALUES ) ' 

( RHM  ) 

6/,'  6)  TUBE  RADIUS/SAMPLE  RADIUS  (B/A) 

7/  ,  '  7)  X  VALUE  ( KDA ) ' 

8/,'  8)  MIN. MAX, STEP  (X  VALUES)' 

9/, '  9)  STORE  DATA' 

9/,'  10)  EXIT  FROM  PROGRAM') 

0121 

N=N-1 

0122 

READ  (5,640)  KO 

0123 

640 

FORMAT  (14) 

0124 

GO  TO  (110.130.150,180,230,270,210.660,650,690)  KO 

0125 

650 

CALL  FILE(N.A.B) 

0126 

GO  TO  620 

0127 

660 

WRITE  (5,670) 

0128 

670 

format;/, '*enter  x  values:  start(known) ,end»step(+/ 

0129 

READ  (5.680)  AKDA , AKM , DKDA 

0130 

680 

FORMAT  (3F15.0) 

0131 

GO  TO  230 

0132 

490 

CONTINUE 

0133 

END 
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SUBROUTINES 


SUBROUTINE: 

DET 

c 

c 

SUBROUTINE  DET 

c 

DECEMBER!  1980 

c 

EDITED  BY  TINA  RUGGIERO 

c 

OOOl 

SUBROUTINE  DET(N.AfD) 

c 

c 

THIS  SUBROUTINE  USES  THE  METHOD  OF  GAUSSIAN  ELIMINATION 

c 

TO  CALCULATE  THE  DETERMINANT 

c 

0002 

DOUBLE  PRECISION  A < 5 . 5 > . B < 5 > » X ( 5 > » D > SUM , C 

0003 

DIMENSION  JPRM ( 5 ) 

0004 

D=1.0D  00 

0005 

DO  13  I *  1 ! N 

0006 

X  ( I ) =0 . OD  00 

0007 

13 

JPRM ( I ) = I 

c 

c 

FIND  THE  ELEMENT  OF  MAXI MUN  ABSOLUTE  VALUE 

c 

0008 

DO  1  K  =  1  i  N 

000? 

C=A(K!K> 

0010 

1 1  =  K 

0011 

JJ=K 

0012 

DO  2  J  =  K  f  N 

0013 

DO  2  I=KiN 

0014 

IF  <DABS(C)-DABS(A(If J) > )3f2f2 

0015 

3 

C=A( I , J) 

0016 

1 1  =  I 

0017 

JJ  =  J 

0018 

n 

CONTINUE 

0019 

D  =  B*C 

0020 

IF ( DABS(D) >20f20f30 

0021 

20 

URITE(SfIOO) 

0022 

100 

FORMAT  < '  MATRIX  A(NfN)IS  SINGULAR') 

0023 

D  =  0 . 0 

0024 

RETURN 

0025 

30 

B  ( 1 1 >  =6  < 1 1 ) /C 

0026 

KPO-K+1 

C 

c 

DIVIDE  EACH  ELEMENT  OF  THE  II  TH  ROW  BY  C 

c 

0027 

DO  4  J=K.N 

0028 

4 

A< I I f J)*A( II , J)/C 

002? 

A  (  I  I • J J  > ” 1 . OD  00 

0030 

IF  <  I  I . EQ.K )G0  TO  60 

c 

c 


SWITCH  THE  KTH  ROW  AND  THE  II  TH  ROW 


0031 

032 

0033 

0034 

0035 

0036 

0037 

0038 

0039 


0040 

0041 

0042 


0043 

0044 

0045 

0046 

0047 

0048 

0049 

0050 

0051 

0052 

0053 

0054 

0055 


C 

D=-D 

DO  5  J=K»N 
C=A ( K  *  J ) 

A(K»J)=A(II»J) 

5  A(II,J)=C 
C  =  B<  K ) 

B(K )=B( II > 

B ( I I ) =C 

60  IF  <  J J . EQ . K ) GO  TO  70 

C 

C  STORE  THE  LOCATION  OF  THE  MAX  PIVOT 

C 

I I=JPRM  <  JJ ) 

JPRM  <  J J ) =JPRM(K) 

JPRM  <  K )  =  1 1 

C 

C  SWITCH  THE  KTH  AND  THE  JJTH  COLUMNS 

C 

D  =  -D 

DO  6  1=1 »N 
C=A ( 1 1 K  > 

A(I»K)=A(I>JJ) 

6  A(I»JJ)=C 

70  IF ( K . EQ . N ) GO  TO  1 

DO  7  I=KPO»N 
DO  8  J=KPO » N 

8  A(I»J)=A(I.J)~A<I»K)*A(K»J) 

B<I)=B<I)-A<I»K>*B(K> 

7  A<IfK)=O.OD  00 

1  CONTINUE 

END 
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SUBROUTINE:  DI3 


0001  SUBROUTINE  DI 3 ( BAH  , ANS  ) 


C 

C  THIS  SUBROUTINE  COMPUTES  2  OPTIONS. 

C  OPTION  1  IS  SOLID-FLUID-RIGID  BOUNDARIES 

C  OPTION  2  IS  SOLID-FLUID-FREE  BOUNDARIES 

C 


0002 

COMMON  AKA. AKDA.BAMS. ACC , COCD , BDA » RHM » KB 

0003 

IF  (BAM)  31.32.32 

0004 

31 

AKA*0 . 

0005 

GO  TO  33 

0006 

32 

AK4=AKDA/BAM 

0007 

33 

AI,SA  =  AKDA/BAMS 

0008 

Q2A=AKA**2-AKDA**2 

0009 

S2A=AKA**2-AKSA**2 

0010 

IF  (RHM)  65.65,40 

0011 

40 

AKFA»AKBA/COCD 

0012 

0F2A=AKA**2-AKFA#*2 

0013 

IF  ( QF2A )  61.60.60 

0014 

60 

QFA=SQRT ( QF2A ) 

0015 

QFB=BDA*QFA 

C 

COMPUTES  BESSEL  FUNCTIONS 

0016 

CALL  BESK(QFA.O.rOO.IER) 

0017 

CALL  B£SK(QFA.1.Y1Q»IER) 

0018 

CALL  BESK(QFB.O.YOOB.IER) 

0019 

CALL  BESK(QFB.l.YlOB.IER) 

0020 

CALL  IO(QFA.BFOQ) 

0021 

CALL  IO(QFB.BFOQB) 

0022 

CALL  INUE(QFAil.BFOQ.BFlO) 

0023 

CALL  INUE(QFB.l.BFOQB.BlOB) 

0024 

QFA»-QFA 

0025 

QFB=-QFB 

0026 

YOQ=-YOQ 

0027 

YOQB*-YOQB 

0028 

GO  TO  65 

0029 

61 

C1FA  =  SQRT  (-QF2A) 

0030 

QFB=BDA*QFA 

0031 

CALL  BESJ(QFA,O.BFOQ,ACC.IER) 

0032 

CALL  BESJ(QFA, 1.BF1Q, ACC. IER) 

0033 

CALL  BESJfOFB.O.BFOQB.ACC. IER) 

0034 

CALL  BESJ(QFB. 1.B1GB.ACC. IER) 

0035 

CALL  BESY(QFA,O.YOO.IER> 

0036 

CALL  BESY(GFA.l.YlO.IER) 

0037 

CALL  BESY(QFB,O.YOQB, IER) 

0038 

CALL  BESY(QFB.l.YlOB.IER) 

0039 

65 

IF(Q2A)8l .80.80 

0040 

80 

0 ASSORT ( Q2 A  > 

004  1 

CALL  IO(GA.BOQ) 

0042 

CALL  INUE(QA.l.BOQ.BlO) 

004  3 

QA= -OA 

0044 

GO  TO  91 

0045 

81 

GA=SGRT(-Q2A> 

0046 

CALL  BESJ(QA,0,BOQ> ACC, IER)  . 

D 

WRITE  (5.85)  IER 

0047 

85 

FORMAT  ( '  IER  =1  ' . I2> 

0048 

CALL  BESJ(QA> 1 »BIQ. ACC. IER) 

0 

WRITE  < 5 * 85 )  IER 

0049 

91 

IF  (S2A)  71 >70  »  70 

0050 

70 

SA*SQRT <S2A) 

0051 

CALL  10 ( SA  > BOS ) 

0052 

CALL  INUE<SA>1>B0S>B1S) 

0053 

00  TO  101 

0054 

71 

SA-SQRT<-S2A) 

0055 

D 

CALL  B£SJ(SA>0>B0S>ACC>IER) 

WRITE  ( 5 f 85 >  IER 

0056 

D 

CALL  BESJ(SAtliBlSiACCrlER) 

WRITE  <  5  >  85 )  IER 

0057 

101 

DU>B0Q*<2.*AKA**2-AKSA**2>+B1Q*2.*QA 

0058 

D12=2.*AKA*(B1S-SA*B0S) 

0059 

D21=2 ,*AKA*QA*B1Q 

0060 

D22=B1S*(2.*AKA**2-AKSA**2> 

0061 

IF  < RHh >  21 >21 r 22 

0062 

21 

ANS=D11*D22-B12*D21 

0063 

00  TO  11 

0064 

D13«AKSA*#2*RHH*BF0Q 

0065 

B14=AKSA**2*RHN*Y0Q 

0066 

D31=QA*B1Q 

0067 

D32-AKA4B1S 

0068 

D33»-QFA*BF1Q 

0069 

D34--QFA*Y1Q 

0070 

IF  (KB-1)  75 >  75  >  76 

0071 

76 

D43»BF0QB 

0072 

D44-Y0QB 

0073 

GO  TO  78 

0074 

75 

D43«<JFB*B1QB 

0075 

D44=QFB*Y1QB 

0076 

78 

P12  =  D11*D22-D12*I>21 

0077 

P34=D33*D44-D34*B43 

0078 

P23=D21*D32-D22*D31 

0079 

P14=D13*044-D14*D43 

0080 

ANS=P12*P34+P23*P14 

0081 

11 

RETURN 

0082 

END 

SUBROUTINE :  CALIB 


OOOl 


0002 

0003 

0004 

0005 

0006 

0007 

0008 

0009 

0010 

0011 

0012 

0013 

0014 

0015 

0016 

0017 

0018 

0019 

0020 

0021 

0022 

0023 

0024 

0025 

0026 

0027 

0028 

0029 

0030 

0031 

0032 

0033 

0034 

0035 

0036 

0037 

0038 

0039 

0040 

0041 

0042 

0043 

0044 

0045 

0046 

0047 

0048 

0049 


c  SUBROUTINE  CALIB ( BAM » ANS ) 

C  THIS  SUBROUTINE  CORRESPONDS  TO  THE  OPTIONS  3»4r  AND  5. 

C  OPTION  3  IS  UHERE  FLUID  SOLID  AND  FREE  BOUNDERIES  IS  ANALYZED 

C  OPTION  4  IS  UHERE  FLUID-S0LID-R10ID  BOUNDERIES ( CEMENTED) 

C  OPTION  5  IS  UHERE  FLUID-SQLID-RIGIDCNOT  CEMENTED) 

C  IS  ANALYZED. 

DOUBLE  PRECISION  AA(StS) 

COMMON  AKA  »  AKDA r  B AMS  >  ACC  >  COCD » BDA  * RHM , KB » AA 
IF <  BAM ) 31 1 32 >  32 

31  AKA=0 

GO  TO  33 

32  AKA=AKDA/BAM 

33  AKSA=AKDA/BAhS 
Q2A=AKA**2-AKDA**2 
S2A»AKA«*2-AKSA**2 
AKF A* AKDA/ COCD 
QF2A=AKA**2-AKFA*#2 
AKB=AKA*BDA 
AKSB=AKSA*BDA 
Q2B«Q2A*BDA*BDA 
S2B=S2A*BDA*BDA 

IF ( QF2A) 61 » 60»  60 

60  QFA»S0RT<QF2A> 

C  BESSEL  FUNCTIONS  ARE  COMPUTED  HERE 

CALL  IO(QFA.BFOG) 

CALL  INUE(QFA*1» BF0Q.BF1Q) 

QFA--QFA 
GO  TO  65 

61  QFA=SQRT (-QF2A) 

CALL  BESJ<QFA»0*BFOQ»ACC»IER> 

CALL  BESJ(QFA>lrBFlQ>ACC>IER) 

65  IF  (Q2A) 81 1 80  *  80 

80  QA*SQRT ( Q2A  > 

QB=S0RT(02B) 

CALL  IO(QA.BOQ) 

CALL  I0( OB  >  BBOQ ) 

CALL  INUE<QA.1»B0Q»B1Q> 

CALL  INUE(QBtlr  BBOQ  r  BB10 ) 

CALL  BESMQAfO.BOArlER) 

CALL  BESK(QBtO>BOBtIER) 

CALL  BESK<QArl>BlArIER) 

CALL  8ESN<QBrl»BlB.IER> 

QA=-QA 
OB=-QB 
BOA=-BOA 
BOB=-BOB 
GO  TO  91 

81  QA=S0RT(-Q2A) 

QB=SQRT< -Q2B) 

CALL  BESJ(QA»0rB0Q>ACC» IER) 

CALL  BESJ<QB.0.BB0Q. ACC. IER) 

CALL  BESJ(QArlrBlQ.ACC>IER) 

CALL  BESJ(QB»l»BB10f ACC.IER) 

CALL  BESY<QA.OfBOA» IER) 

CALL  BE3Y<Q8.0fB0BfIER> 


0050 

CALL  BESY < QA *  1  * B1 A* IER > 

0051 

CALL  BESY(QB* 1 *B1B* IER) 

0052 

91 

IF<S2A>71*70*70 

0053 

70 

SA«SQRT <S2A) 

0054 

SB*SQRT(S2B> 

0055 

CALL  IO<SA*BOS) 

0056 

CALL  I0(8B*BB0S) 

0057 

CALL  INUE<  SA* 1  *  BOS  >  B IS ) 

0058 

CALL  I HUE (SBfIfBBOSfBBIS) 

0059 

CALL  BESK(SA*0*BSA0*ZER> 

0060 

CALL  BE8K(SB*0*B8BO*IER> 

0061 

CALL  BESK(SA*liB8Al*IER> 

0062 

75 

CALL  BESK(8B*1*B8B1*!ER) 

0063 

77 

BSA1— BSA1 

0064 

BSB1  — BSB1 

006  5 

GOTO  101 

0066 

71 

SA“SQRT  < -S2A  > 

0067 

SB*SQRT ( -S2B ) 

0068 

CALL  BESJ(SA*0*B08i ACC > IER > 

0069 

CALL  BESJ(SB*0* BBOS » ACC  * IER ) 

0070 

CALL  BESJ<SA*liBlS*ACCiIER> 

0071 

CALL  BE8J(SB* 1»BB1S*ACC*IER> 

0072 

CALL  BE8Y(SA*0>BSA0* IER) 

0073 

CALL  BESY(SBfO*BSBO*IER> 

0074 

CALL  BESYtSAil* B8A1 * IER ) 

0075 

85 

CALL  BESY(SB*1*BSB1*IER) 

0076 

87 

CONTINUE 

C 

THE  DETERMINANT  IS  BEING  FILLED 

0077 

101 

AA<1>1>-B0Q*<2.*AKA**2-AKSA**2>+2.*(QA*B1Q> 

0078 

AA( 1f2>»(2.*AKA**2-AKSA**2)*B0A+2.*QA*B1A 

0079 

AA  < 1  *  3) *2 • *AKA*< B1S-9A4B0S ) 

0080 

AA ( 1  *  4 ) -2 . * AKA* ( BSA1 -S A*BSA0 ) *6A 

0081 

AA( 1 .5>-RHM*AKSA*AKSA*BF0Q 

0082 

AA<  2  > 1 ) “2 . *AKA*QA*B10 

0083 

AA ( 2  f  2 ) "2 . *AKA*QA*B1A 

0084 

AA  <  2 . 3 ) - <  2 . *AKA*AK A- AKS A* AKSA >  *B1S 

0085 

AA <  2  *  4 ) - <  2 . *AK A* AK A-AKSA*AKSA ) *BS A 1 *SA 

0086 

AA ( 2  f  5 ) *0 . 0 

0087 

IF < KB  .EG.  3) GO  TO  120 

0088 

AA ( 3  > 1 ) «BB 1Q*QB 

0089 

AA(3.2)*B1B*QB 

0090 

AA(3f3)-AKB*BB1S 

0091 

AA<3f4)»AKB*BSB1*SA 

0092 

GO  TO  150 

0093 

120 

AA(3f1)sBB0Q*<2f*AKB**2-AKSB**2)42>*QB*BB1Q 

0094 

AA(3f2)«(2.*AKB*AKB-AKSB*AKSB)*B0B+2.*QB*B1B 

0095 

AA(3>3)-2.*AKB*(BB1S-SB*BB0S) 

0096 

AA(3,4)»2.*AKB*(BSB1-SB*BSB0)*SA 

0097 

150 

AA(3f5)»0.0 

0098 

IF ( KB  .NE.  4 ) GO  TO  160 

0099 

AA(4*1)»AKB*B80Q 

0100 

AA(4f2)»AKB*B0B 

0101 

AA(4f3)*-SB*BB0S 

0102 

AA(4f4)*-SB*BSB0*SA 

0103 

GO  TO  170 

0104 

160 

AA(4*1)*2.*AKB*QB*BB1Q 

0105 

AA(4f2)»2.*AKB*QB*B1B 

0106 

AA<  4*3)«(2.*AKB*AKB-AK8B*AKSB>*BB1S 

0107 

AA<4f4)»(2.*AKB*AKB-AKSB*AKSB)*BSB1*SA 

0108 

170 

AA( 4  f  5) “0 . 0 

0109 

AA  <  5 • 1 ) »Q AtBlQ 

0110 

AA  <  5  >  2 ) «QA*B  1 A 

0111 

AA(5f3)»AKA*B1S 

0112 

AA(5f4)“AKA*BSA1*8A 

0113 

AA<5f5)«-QFA*BF10 

0114 

IF  < RHM  .EQ.  0.0)AA(5f5)«1.0 

0115 

CALL  DET ( 5 • AA • AN8  > 

0116 

RETURN 

0117 

END 
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SUBROUTINE:  FILE 


0001 

SUBROUTINE  FILE <N » FREQ. RATIO ) 

0002 

BYTE  FIL£X<32) 

0003 

REAL  FREQ< 1000) »RATIO< 1000) 

0004 

WRITE (5 1 10 ) 

0003 

10 

FORMAT (/» 'BFILE  NAME  FOR  X.Y  DATA?  ') 

0006 

READ(5f 20)LENtFILEX 

0007 

20 

FORMAT  (Q.32A1 ) 

0008 

FILEX (LEN  4  1 > -0 

0009 

OPEN< UNIT-2. NAME-FILEX. TYPE- 'NEW' ) 

0010 

DO  40  J-l.N 

0011 

WRITE<2»30)FREGKJ) .RATIQC J) 

0012 

30 

FORMAT (2E20. 10) 

0013 

40 

CONTINUE 

0014 

CLOSE ( UNIT-2 ) 

0015 

CONTINUE 

0016 

WRITE  <  5 . 43 ) 

0017 

45 

FORMAT <  / » '  FILE  IS  STORED  IN  1  FILE  WITH  FORMAT  (2E20.0) 

0018 

RETURN 

0019 

END 
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TEST  METHODS  AND  RESULTS 


RDM  UGUIDE  <CR> 

ENTER  *  CORRESPONDING  TO  OPTION  DESIRED 

1)  SAMPLE  IN  FLUID  TUBE  UITH  RIGID  BOUNDARY 

2)  SAMPLE  IN  FLUID  TUBE  UITH  FREE  BOUNDARY 

3)  FLUID  FILLED  ELASTIC  TUBE.  OR  EMPTY  TUBE 

4)  FLUID  FILLED  COATED  RIGID  TUBE  (CEMENTED) 

5)  FLUID  FILLED  COATED  RIGID  TUBE  (NOT  CEMENTED) 

6)  FLUID  TUBE  UITH  RIGID  BOUNDARY 

7)  FLUID  TUBE  UITH  FREE  BOUNDARY 
3. <CR> 

ENTER  RATIO  OF  SHEAR  TO  BULK  MODULUS  (G/K)  .  279  <CR> 

ENTER  RATIO  OF  FLUID  TO  SOLID  SPEED  (CO/CD)  .315  <CR> 

ENTER  RATIO  OF  FLUID  TO  SOLID  DENSITY 
(FOR  EMPTY  TUBE  SET  RHM=0, VACUUM  TUBE  SET  RHM*NEG  VALUE)  0  <CR> 

SHE AR=  .431E+00  YOUNG=  .747E+00  KAYLEIGH=  .423E+00 

ENTER  X  VALUE  (SET  KDA  NEGATIVE  FOR  ASYMPTOTES)  -1  «CR> 

ENTER  MIN. MAX. STEP  (HORIZONTAL  SCAN)  ■  1  .  4  .  .  .05  <CR> 

ENTER  RATIO  OF  TUBE  TO  SAMPLE  RADIUS  (B/A)  1.134  <CR> 

ASYMPTOTE  AT  X  =  0.73625  4/-  0.230E-03 

NO  MORE  ASYMPTOTES 

ENTER  *  CORRESPONDING  TO  CHANGE 

1)  OPTION 

2)  SHEAR  MODULUS/BULK  MODULUS  (G/K) 

3)  SPEED  IN  FLUID/SPEED  IN  SOLID  (CO/CD) 

4)  FLUID  DENSITY/SOLID  DENSITY  ( RHM  > 

5)  MIN  MAX.  STEP  <  Y  VALUES) 

6)  TUBE  RADIUS/SAMPLE  RADIUS  'B/A> 

7)  X  VALUE  (KDA) 

3)  MIN. MAX. STEP  (X  VALUES) 

9)  STORE  DATA 

10)  EXIT  FROM  PROGRAM 
8  <CR> 

ENTER  X  VALUES'.  START  (KNOUN)  .  END  .  STEP  (  +  /  -  /  ■  t  •  2 .5  .  .5  <CR> 

ENTER  MIN. MAX, STEF  (VERTICAL  SCAN)  .  1  .  1  .  5  ,  .  1 
X  =  0.10000  Y  =  0.74605  ERROR  =  0.500E-04 

X  *  0.60000  Y  =  0 . 62565  ERROR  =  0.500E-0-1 

X  -  1.10000  v  =  0.19825  ERROR  *  0.500E-04 

<  *  1.60000  Y  =  0.21645  ERROR  =  0.300E-04 

Y.  *  2.10000  Y  =  0.23655  ERROR  =  0.500E-04 

ENTER  *  CORRESPONDING  TO  CHANGE 

1)  OPTION 

2)  SHEAR  MODULUS/BULK  MODULUS  (G/K) 

3)  SPEED  IN  FLUID/SPEED  IN  SOLID  (CO/CD) 

4)  FLUID  DENSITY/SOLID  DENSITY  (RHM) 

5)  MIN  MAX.  STEP  (Y  VALUES) 

6)  TUBE  RADIUS/SAMPLE  RADIUS  (B/A) 

7)  X  VALUE  (KDA) 

3)  MIN. MAX. STEP  (X  VALUES) 

9)  STORE  DATA 

10)  EXIT  FROM  PROGRAM 

a  <CR> 
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♦FILE  NAME  FOR  X.Y  DATA: 

DATA  <CR> 

ENTER  ♦  CORRESPONDING  TO  CHANGE 

1)  OPTION 

2)  SHEAR  MODULUS/BULK  MODULUS  (G/K) 

3)  SPEED  IN  FLUID/SPEED  IN  SOLID  (CO/CD) 

4)  FLUID  DENSITY/SOLID  DENSITY  ( RHM ) 

.5)  MIN  MAX  >  STEP  (Y  VALUES) 

6)  TUBE  RADIUS/SAMPLE  RADIUS  (B/A) 

7)  X  VALUE  <  KDA  > 

8)  MIN  » MAX  *  STEP  (X  VALUES) 

9)  STORE  DATA 

10)  EXIT  FROM  PROGRAM 
10  <CR> 


l)  All  underlined  portions  are  user  supplied. 

2)  <CR>  indicates  'RETURN' 


5 

i 
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DESCRIPTION  or  PROGRAM  IMPED 


The  program  IMPED  computes  values  £or  Che  complex  dimensionless 
propagation  speed  in  Che  solid  sample  (assumed  infinite  in  length)  inside  the 
impedance  tube.  Input  to  the  program  are  the  complex  values  for  the  ratio  of 
shear  modulus  to  bulk  modulus  G/K,  and  the  speed  in  the  fluid  divided  by  the 
dilatational  wave  speed  in  the  solid,  c  /cd.  Thus,  it  may  be  considered  as 
the  complex  counterpart  to  program  WGUIDE,  Option  1.  The  root  finding 
technique  is  identical  to  that  used  in  program  KTUBE.* 

Subroutines  for  Bessel  functions  were  written  by  direct  series  expansion 
since  the  real  part  of  the  argument  is  small  enough  that  a  small  number  of 
terms  suffices.  The  names  of  the  subroutines  are  CBJO,  CBJ1,  CBYO,  and  CBY1, 
where  the  last  two  letters  of  each  name  indicate  kind  and  order  of  the  Bessel 
function. 

The  output  produced  by  the  program  is  in  the  form  of  real  and  imaginary 
parts  of  the  relative  wave  speed  S  *  c/cd,  each  with  its  own  error  limit, 
corresponding  to  the  input  error  limits  set  separately  for  the  two  terms.  The 
attenuation  constant  may  be  readily  calculated  by  a  ■  -k^B  /(B^  +  B-2) 
where  a  is  the  attenuation  per  unit  of  length  and  Bj  andp^  are  the  real  and 
imaginary  parts  of  B  which  is  equal  to  c/cd.  If  the  propagation  of  the  wave 
in  water  is  assuawd  to  be  without  loss  (cQ  real),  the  input  value  of  AKDA 
should  have  the  same  ratio  of  imaginary  to  real  part  as  COCD. 


Important  Variables  in  IMPED 

Input  Variables 

AKDA  (k.a),  Dimensionless  wave  number  of  dilatational  waves, 

•  <fua)/cd  (real  &  complex) 

AMI  #  of  iterations 

AMR  #  of  divisions  of  imaginary  part 

BAM  (c/cd)  relative  wave  speed;  enter  real  part  of  c/cd  (ZNOT) 

BDA  (b/a)  b  *  radius  of  tube 

a  ■  radius  of  cylindrical  sample 

COCD  (cQ/cd)  ”  wave  sP®ed  in  fluid  divided  by  dilatational  wave 

speed  in  solid  (real  and  complex) 

DELZ  (used  in  computation  of  steps ize) 

GOK  (G/K)  (real  and  complex)  G  ■  shear  modulus,  K  *  bulk  modulus 


*  P.S.  IXibbelday,  "Complex  Root-Finding  Program  with  Application  to  the 
Dispersion  Relation  of  Waves  Propagating  in  a  Fluid  Plate,"  NRL  Memorandum 
Report  4559,  November  1981. 
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RHM 

T0L1 

TOL2 

Computat i 
AI 
AKA 

AKDA 

AKFA 

AKSA 

ANR 

AMGA 

BAMS 

DEL 

DELH, 

DELP 

DELZ 

DFDZ 

DS(Z) 

DZ 

FI 

F2 

F3 

F4 


(p0/pg)  ratio  of  density  of  fluid  to  density  of  solid;  if 
RHM  <  0  solid  is  in  vacuum 

Tolerance  of  real  part 

Tolerance  of  imaginary  part 


>nal  Variables 


i,  imaginary  unit 

(ka),  dimensionless  wave  number,  * — 
a  *  radius  of  cylinder 
(i)  ”  angular  frequency 
c  =  propagation  speed 

(k^a)  Dimensionless  wave  number  of  dilatational  waves, 

=  ua/c^,  c^  =  dilatational  wave  speed 

k  a  •  Dimensionless  wave  number  in  fluid 
o 

(kga),  Dimensionless  wave  number  of  shear  waves  kga  *  ma/cs 

#  of  divisions  of  imaginary  part 

Imaginary  part  of  wave  speed  8 

(cg/cd)»  relative  wave  speed  of  shear  waves 

Real  stepsize 

DELV  Step  sizes  for  movement  of  test  pairs,  horizontally  and 
vertically. 

Intermediate  symbol  of  stepsize  DEL. 

Az,  used  in  computation  of  approximation  to  3f/3z 
Approximation  for  df/dz. 

Function  subprogram 

Complex  stepsize 

Function  evaluated  at  Zl:  f(z^) 

Function  evaluated  at  Z2:  f(z2) 

Function  evaluated  at  Z3:  f(zj) 

Function  evaluated  at  Z4 :  f(z^) 
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FV1,  FV2 

FH1 ,  FH2 

Flags  to  Indicate  direction  of  motion  of  test  pairs 

FSTEP 

Factor  to  adjust  step  size  DEL 

IOP 

Return  option  number 

■  1  Program  returns  to  entering  COCD,  GOK,  AKDA 
=  2  Program  returns  to  entering  Real  Seed,  #  of 

Iterations,  DELZ,  BDA 

*  3  Program  returns  to  entering  AWR,  TOL  (Real),  TOL 
(Imag.),  RHM 

*  4  Program  exits 

KAM 

-  AMI 

KR,  KL 

KUP,  KDOWN 

Counters  in  loops  to  check  number  of  iterations 

NR 

-  ANR 

NRV 

Counter  in  advancing  the  parameter  RHO 

Q2A 

(qa)2  -  (ka)2  -  <kda)2 

QF2A 

(qQa)2  ■  (ka)2  -  (kQa ) 2 

QF2B 

(qQb)2  -  (kb)2  -  (kea)2 

REGA 

Real  part  of  relative  wave  speed  3 

RHM 

Nominal  density  of  fluid  loading  the  plate,  divided 
density  of  plate  material* 

by 

RHO 

Stepwise  varied  value  of  relative  density,  varying 
zero  to  RHM 

from 

SIGN 

Expression,  the  algebraic  sign  of  which  determines 
of  root  relative  to  test  pair. 

location 

TOL 

Limit  for  relative  error  in  the  root. 

Z 

Variable  for  root  used  in  calling  subroutines 

ZI,  Z2 

Vertical  test  pair 

Z3,  Z4 

Horizontal  test  pair 

ZNOT 

Value  for  the  real  root  of  the  dispersion  relation 
real  root  finder,  serving  as  the  ..eed  for  starting 
program. 

from  a 
the 

ZVAR 

Intermediate  value  of  root,  serving  as  seed  for  the 
step,  in  parameter  RHO 

next 
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ZVAR 


Intermediate  value  of  root,  serving  as  seed  for  the  next 
step,  in  parameter  RHO 


Output  Variables 

AMGA  Imaginary  part  of  root 

ZEGA  Real  part  of  root 

Hajor  Computation  Block  in  IMPED 

Refer  to  major  computation  blocks  in  R0819 
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COMPLEX  ROOTFINDING  ROUTINE 
WRITTEN  BY  TINA  RUGGIERO 
DECEMBER  17.1981 


THIS  COMPLEX  ROUTINE  INTRODUCED  BY  DUBBELDAY 
ALLQUS  THE  USER  TO  SOLVE  FOR  BOTH  THE  REAL  AND  IMAGINARY 
PARTS  OF  A  COMPLEX  ROOT.  GIVEN  A  REAL  ROOT  IN  THE 
VICINITY  AS  A  SEED. 


SUBROUTINES  USED! 

CDI3  (CONTAIN  FUNCTION  DESCRIBED  BY  DETERMINANT  OF  A  SET  OF 
BOUNDARY  CONDITIONS) 

CB JO  (BESSEL  FUNCTION  <JO>> 

CBJ1  (BESSEL  FUNCTION  <J1>) 

CBYO  (BESSEL  FUNCTION  (YO)> 

CBY 1  (BESSEL  FUNCTION  (Yl>> 


OOOl 

COMPLEX  Z.DSP.AI.DFDZ.DZ.ZVAR.Z1.Z2.Z3.Z4.F1.F2.F3.F4. 

1  GOK  r  DSP2 • CDI3 .CMF'X . AKDA. AKA >  BAMS » COCD . HGOK  > HAKDA >  HCOCD 

0002 

P 

COMMON  AKA. AKDA, COCD. BDA.RHM. GOK 

L 

C 

INPUT  SECTION 

0003 

5 

WRITE ( 5  » 1 0 ) 

0004 

10 

FORMAT ( / ' *ENT£R  COCD ( R » I ) , GOK ( R » I > , AKDA ( R , I )  ') 

0005 

READ (5, 20) HCOCD, HGOK. HAKDA 

0006 

20 

FORMAT (6F10.0) 

0007 

25 

UR  I TE  <  5 , 30  ) 

0008 

30 

FORMAT (/'*ENTER  REAL  SEED,*  OF  ITER . , DELZ , BDA  ') 

0009 

READ(5»40)ZN0T.AMI,DELZ,BDA 

0010 

40 

FORMAT (4F10.0) 

0011 

K AM1 AM I 

0012 

BAMS=CSQRT(G0K/(l.+4.*G0K/3. ) ) 

0013 

45 

URITECS, 50 ) 

jO  1  4 

50 

FORMAT (/'CENTER  ANR.TOL  ( REAL ) , TOL ( IMAG > ,RHM  ') 

0015 

READ (5. 60) ANR.T0L2.T0L1 »RHM 

0016 

60 

F  ORMAT ( 4F 1 0 • 0 ) 

0017 

AI=CMPLX(0. ,1.0) 

0018 

GOK=HGOK 

0019 

COCD-HCOCD 

0020 

AKDA=HCOCD 

0021 

RGOK=REAL( GOK) 

0022 

R AKDA  =  RE AL ( AKDA  > 

0023 

YGOK=AIMAG(GOK) 

0024 

YAKDA=AIMAG(ANDA) 

0025 

RCOCD=REAL(COCD) 

0026 

YCOCD=AIMAG(COCD> 

0027 

NR= I NT ( ANR  ) 

0028 

NR  V  =  0 

0029 

Z V AR  =  ZNO  T 

0030 

GOK=CMPLX(RGOK.O.O) 

0031 

AKDA =CMPLX(R AKDA, 0.0) 

0032 

C0CD=CMPLX(RC0CD,0.0) 

0033 

Z=ZVAR+DELZ 

0034 

DFDZ=CDI3(Z)/DELZ 
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GOK»GOK+AI*YGOK/ANR 

AKDA*AKDA+AI*YAKDA/ANR 

COCD=COCD+AI*YCOCD/ANR 

NRV-NRV+1 

DZ=-CDI3(ZVAR)/DFDZ 
X»ABS( REAL ( DZ ) ) 

Y«ABS( AlhAG(DZ)  ) 

DELI  AND  DEL2  ARE  HORIZONTAL  AND  VERTICAL  STEP  SIZES 

DEL1=Y 

DEL2=X 

WRITE(5. 16) DELI >  DEL2 

FORMAT  ( / »  '  DELI*  '.£8.2.'  DEL2*  ',EB.2> 

WRITE ( 5  » 57 ) 

TYPE  *.' - 

DELP=DEL1 

DELQ-DEL2 

FIND  DIRECTION  OF  INITIAL  VERTICAL  PAIR 

DEL2=DELP 

DEL1*DELQ 

Z1=ZVAR-AI*DEL1 

Z2=ZVAR+AI*DEL1 

F 1=CDI3  <  Z1 ) 

F2=CDI3(Z2> 

KR  =  0 
KL  =  0 
KUP  =  0 
KD0WN=0 

THE  RESULT  OF  SIGN  IS  TO  DETERMINE  DIRECTION  OF  PAIR 
SIGN=AIMAG(F1 > *RE AL ( F2 > -REAL < F 1 > *A IMAG ( F2 > 

IF  (SIGN)  71.72.73 

FH1=-1 
GO  TO  74 

TYPE  # . ' LEFT  RIGHT  SIGN  ' 

GO  TO  100 
FH 1  =  1 

DELH=FH1*DEL1 

LOOP  TO  FIND  VERTICAL  LINES  BRACKETING  ROOT 

Z1=Z1+DELH 
Z2=Z2+DELH 
FI =CDI3  <  Z1 ) 

F2=CDI3  <  Z2 ) 

SIGN  IS  USED  TO  DETERMINE  DIRECTION  PAIR  IS  MOVED 
SIGN*AIMAG(F1 > *RE AL ( F2 ) - A I  MAG ( F2 > *RE AL ( F 1  ) 

IF  (SIGN)  81.72.83 

FH2=-1 

KL=NL+ 1 

IF (KAM-KL)9?.99»S4 
TYPE  *.  'EXIT  LEFT' 

GO  TO  100 


0081 

IF(KAH-KR)98>98.84 

0082 

98 

TYPE  *. 'EXIT  RIGHT ' 

0083 

GO  TO  100 

0084 

84 

IF<FH1*FH2>92.75.75 

008S 

C 

FIND  HORIZONTAL  TEST  PAIR  DIRECTION 

92 

Z4=(Zl+Z2)/2.+DEL2 

0086 

KR«0 

0087 

KL=0 

0088 

Z3=(Zl+Z2)/2. -DEL2 

0089 

F3=CDI 3 ( Z3 ) 

0090 

F4=CDI3(Z4) 

0091 

S1GN»AIMAG<F3)*REAL<F4)-AIHAG<F4>*REAL(F3> 

0092 

IF(SIGN) 101 . 102. 103 

0093 

101 

FU1  =  1 . 

0094 

GOTO  104 

0095 

102 

TYPE  * « ' UP- DOWN  SIGN=  ZERO' 

0096 

GOTO  100 

0097 

103 

FU1=-1 

0098 

104 

D£LV=FY1*DEL2 

C 

LOOP  TO  FIND  HORIZONTAL  LINES  BRACKETING  ROOT 

0099 

105 

Z4=Z4+AI*DELU 

0100 

Z3=Z3+AI 4DELV 

0101 

F3=CDI3(Z3> 

0102 

F4=CD 1 3  <  Z 4 ) 

0103 

SIGN»AIMAG(F3)*REAL(F4)-AIMAG(F4)*REAL(F3) 

0104 

IF(SIGN) 111 >  102 i  113 

0105 

111 

FU2  =  1 . 

0106 

KUP=KUP+ 1 

0107 

n 

IF(KAM-KUP>97.97» 114 

0108 

L 

97 

TYPE  *.  'EXIT  UP' 

0109 

GOTO  100 

0110 

113 

FV2*-1  . 

0111 

KD0UN=KD0UN+1 

0112 

IF (KAM-KDOUN 196.96. 114 

0113 

96 

TYPE  *.  'EXIT  DOUN' 

0114 

GO  TO  100 

0115 

114 

IF (F014F02) 122. 105. 105 

0116 

122 

REGA  =  REAL ( Z 1  > 

0117 

AMGA=AIMAG(Z3) 

0118 

KUF=0 

0119 

KD0UN=0 

0120 

133 

AMGA=ABS<AMGA) 

0121 

DREG=ABS(REGA-ZNOT) 

C 

c 

CHECK  IF  REQUIRED  PRECISION  REACHED 

0122 

IF  <  DEL  1 / AMGA-TQL 1  ,LE.  01G0T0  132 

0123 

IF (DEL2/REGA-T0L2  .LE.  OIGO  TO  135 

0124 

DEL2  =  DEL 2/  1 0  . 

0125 

GO  TO  92 

0126 

135 

DEL  1  =  BEL1/ 1 0 . 

0127 

DEL2=DEL2/10. 

0128 

Zl=(Z3+Z4)/2.-AI«DELl 

0129 

Z2=(Z3+Z4)/2.+AI*DELl 

0130 

GOTO  65 

0131 

132 

REGA=REAL ( Z1 ) -DELH/2. 
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F.: 

PI 


b 


ZSZ' 


vt 

•■•I- 


0132 

AMGA=AIMAG(Z3)-0ELU/2. 

0133 

c 

DELH=ABS<D£LH/2.  ) 

0134 

L 

r 

D£LV*ABS(DELV/2.  ) 

u 

c 

OUTPUT  SECTION 

0135 

57 

FORMAT ( / » '  REAL' .  7X.  ' IMAGINARY' ,5X» '1 

0136 

WRITE<5»58)REGA» AMGA . DELH > DELU 

0137 

58 

FORMAT  <E10.4>4X»E10.4»4X»E10.4»6X>E10.4 

C 

C 

CHECK  IF  FINAL  GOK  AND  AKDA  IS  REACHED 

0138 

IF(NR-NRV) 145. 145. 146 

0139 

146 

ZUAR=REGA+AHGA*AI 

0140 

GOK=GOK+AI*YGOK/ANR 

0141 

akda»akda+ai*yakda/anr 

0142 

CQCD*COCD+AI*YCOCD/ANR 

0143 

NRU =NRV+  1 

0144 

GOTO  70 

0145 

145 

CONTINUE 

0146 

TYPE  *.  '  ' 

0147 

100 

TYPE  ENTER  OPTION  1  <1.2. 3. 4) 

0148 

RE AO  <  5 . 160) IOP 

0149 

160 

FORMAT  <121 

0150 

GO  TO  <5.25.45.550) , IOP 

0151 

550 

END 

REAL  ERROR ' *  5X  .  ' IMA6  ERROR' 


« 


4 


« 


0001 

FUNCTION  CDI3 ( BAM > 

0002 

COMMON  AKA.AKDA,COCDrBDA>RHM»GOK 

0003 

COMPLEX  AKAiAKDAiBAMSt COCDt AKSAtQ2A»S2Af AKFA tQF2A>QFArQFB>BF0Q 
1  BFlQrBlQBf TOO. YlQ. Y1QB » QA » BOO. BIO f SA » BOB » B18 » D1 1 > D12» D21 » D22 
1  DI3>D14»D31»B32,D33»D34»D43*D44.P12»P34»P23tP14»CDI3»QP2Af AI 
1  QPF2A » SP2A » QPFA » BAH r D13 r CSQRT >60K 

0004 

BAHS»CSQRT( GOK/ (1.+4.4G0K/3.  >  > 

0005 

AI-CMPLX(O.Oil.O) 

0004 

AKA-AKDA/BAM 

0007 

AKSA-AKDA/BAM8 

0008 

QP2A-AKDAM2-AKAM2 

0009 

SP2A"AKSAM2-AKA«2 

0010 

IF(RHM)65»  45 >  40 

oon 

40 

AKFA-AKDA/COCD 

0012 

QPF2A-AKFAM2-AKAM2 

0013 

IF < REAL (0PF2A)) 45,46.46 

0014 

46 

QPFA-CSQRT(QPF2A) 

0015 

GO  TO  50 

0014 

45 

OPFA"-AI<CSQRT  < -QPF2A ) 

0017 

50 

QFB-BDABQPFA 

0018 

IF ( ABS  <  QFB )  .GT.  10.0)G0T0  1000 

0019 

CALL  CBJO<  QPFA »  BFOQ ) 

0020 

CALL  CBJKQPFAt BF1Q) 

0021 

CALL  CBJl(QFBiBlQB) 

0022 

CALL  CBYO(QPFA.YOQ) 

0023 

CALL  CBYl <QPFA»Y1Q> 

0024 

CALL  CBYHQFB'YIQB) 

0025 

65 

IF ( REAL ( QP2A ) )81»80i80 

0026 

80 

QA"CSQRT ( QP2A  > 

0027 

GO  TO  85 

0028 

81 

QA  — AI<CSQRT<-QP2A> 

0029 

85 

CALL  CBJ0<QA»80Q> 

0030 

CALL  CBJ1 (QAiBlQ) 

0031 

IF<REAL<SP2A)>71»70»70 

0032 

70 

SA*CSQRT (SP2A) 

0033 

GO  TO  90 

0034 

71 

SA=-AI<CSQRT<-SP2A> 

0035 

90 

CALL  CB J0<  SA  » BOS ) 

0036 

CALL  CBJ1 ( SA  r  BIS) 

0037 

D11-B0Q*(2.*AKA**2-AKSA**2)+B1Q*2.*QA 

0038 

D12-2.<AKA<<B1S-SA<B08> 

0039 

D21=2. <AKA<QA<B1Q 

0040 

D22  =  B1S<(2.<AKAM2-AKSA«2) 

0041 

IF  ( RHM ) 21 1 21 r  22 

0042 

21 

CDI3»D11<D22-D12<D21 

0043 

GO  TO  11 

0044 

22 

D13»AKS A<<2<RHM<BF0Q 

0045 

D 14" AKSA<<2<RHM<Y0Q 

0046 

D31*QA<B1Q 

0047 

D32"AKA<B1S 

0048 

D33»-QPFA<BF1Q 

0049 

D34*-QPFA<Y1Q 

0050 

D43"QFB<B1QB 

0051 

D44"QFB<Y1QB 

0052 

P12«D11<D22-D12<D21 

0053 

P34»D33<D44-D34<D43 

0054 

P23  =  D21<I>32-r22<D31 

0055 

P14=D13<D44-D14<D43 

0056 

CDI3*P12<P34+P23<P14 

0057 

11 

RETURN 

0058 

1000 

TYPE  <» 'ARGUMENT  FOR  BESSEL  FUNCTIONSf  TOO  LARGE' 

0059 

2000 

ENO 

0001 


SUBROUTINE  CBJKZ.ANS) 


C 

c 

c 

c 

0002 

0003 

0004 

0005  8 

0006 
0007 
0008 
0009 
0010 
0011 
0012 
0013 
0014 

0015  30 

0016  50 

0017 


BESSEL  FUNCTION  J1 
(CONFUTED  BY  SERIES) 

COMPLEX  Z»FUNC.AN8*CMPLX. FIRST 

CALL  ERRSET <73 »». FALSE. »  *  .FALSE. » 

TOL-. 000001 

ERROR-1.0 

FUNC*Z/2 

ANS-Z/2 

RHO-CABS(Z) 

DO  30  R*2  « 1000  . 

FIRST»-1.0*FUNC*(Z/2.0)**2/(R*R> 
FUNC-FIRST*2.0*R/<2*R-2 .0) 
ANS=ANS+FUNC 

ERROR  =  ERROR*<  RHO/2 . 0)**2/(R*R) 

IF < ERROR/CABS (ANS)-TQL) 50 *50 *30 

CONTINUE 

RETURN 

END 


SUKKOUTINE  CBJO(Z.ANS) 


BESSEL  FUNCTION  JO 
(CONFUTED  BY  SERIES) 

COMPLEX  ZrFUNCtANStCMPLX 
ANS-CMPLX ( 1 • 0  f  0 • 0  > 

TOL-. OOOOOl 
ERROR-1.0 

FUNC-CHPLX ( 1 • 0 i 0 • 0 ) 

RHO-CABS(Z) 

DO  30  R-1.1000 

FUNC- (-1.0) *FUNC* ( Z/2 . 0 ) **2/ ( R*R > 

ANS-ANS+FUNC 

IF  (CABS(ANS)  .EQ.  O.OJGO  TO  30 
ERROR-ERROR* (RH0/2.0) **2/ (R*R) 

IF(£RROR/SQRT ( REAL ( ANS ) **2+AINA6( ANS ) **2 ) -TOL ) 30 * 30 » 30 

CONTINUE 

RETURN 


RUN  IMPED 

ENTER  COCDtR. I > »GOK(R» I > , AKDA(R» I >  . 2358 . . 002358 , . 45 A  1 r  ,  00  454  l  ,  g  ,  ,  QQ2 

ENTER  REAL  SEED.*  OF  I  TER . . OELZ . BD A  97875 . 200  >  .  01  »  1  . 00 1  <CR> 

ENTER  ANR.TOL  ( REAL > » TOL < I  MAG >  .  RHM  3 » . 01 » , QQ1 , . 12805  <CR* 

DELI  *  0.24E-03  DEL2=>  0.16E-05 

ftEAL  IMAGINARY  REAL  ERROR  IMAG  ERROR 


. 9788E+00 

. 9788E+00 

. 9788E+00 

ENTER  OPTION 

4  <CR> 


0 . 23 13E-03 
0.4870E-03 
0.71 83E-03 
*  <1.2. 3. 4) 


0 . 8249E-07 
0 . 8249E-07 
0 . 8249E-07 


0. 1217E-04 
0. 1217E-04 
0. 1217E-04 


1)  All  underlined  portions  are  user  supplied, 

2)  <CR>  indicates  1  RETURN  * 


<CR> 
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APPENDIX  B 

Correction  for  Inperfect  Reflector  in  Impedance  Tube 


The  analysis  of  the  role  of  the  reflector  in  the  impedance  tube  starts 
from  the  observation  that  there  are  traveling  waves  in  both  directions  in  the 
water  of  the  Impedance  tube,  in  the  sample,  and  in  the  reflector,  each  with 
its  own  amplitude  and  phase.  By  equating  stress  and  particle  velocity  at  each 
of  the  interfaces  one  may  infer  the  reflection  coefficient  for  the  following 
cases:  rsa  is  the  reflection  coefficient  without  the  sample  in  place,  rst  is 
the  reflection  coefficient  without  the  sample.  One  would  like  to  relate  these 
measured  quantities  to  the  reflection  coefficient  r  used  in  Eq.  (2),  which 
implies  a  perfect  reflector.  The  length  l  of  the  reflector  is  chosen  such 
that  kj.4  is  close  to  x/2.  The  subscript  r  refers  to  the  reflector.  Define  6 
by  kr&  »  u/2  +  6 ,  6  is  a  small  complex  number.  One  finds  that 

1  -  i(tan  k  d)(pc)  /(pc)  -  ie[l  -  i(tan  k  d)(pc)  /(pc)  ] 

SOS  S  S  O  /‘Di  \ 

fsa  1  +  l(tan  k  d)(pc)  /(pc)  +  ie[l  +  l(tan  kd)(pc)  /(pc)  ]  *  ' 

SOS  so 


r 


st 


e 


-2ik  d  1  -  le 
0  1  +  ie  * 


(B2) 


1  -  i(tan  k  d)(pc)  /(pc) 
_ s _ o  _ s_ 

1  +  i(tan  k  d)(pc)  /(pc) 
s  os 


(B3) 


where  the  parameter  e  is  defined  by  e  ■  (tan  6)(pc)0/(pc)r . 

In  the  present  data  analysis,  the  reflection  coefficient  r  in  Eq.  (2)  is 
replaced  by  rc,  defined  by 

rc  =  exp(-2ikod)raa/rst  .  (B4) 

The  correction  factor  applied  to  rga  according  to  this  equation  may  be 
expressed  in  terms  of  the  small  parameter  e  by 

exp(-2ikQd)/rgt  -  1  +  2ie  .  (B5) 


The  proper  correction  factor  is  r/rga, 
order  by 


which  may  be  expressed  to  the  same 


-  1  +  2iFe 


where  the  factor  F  is  given  by 


F  ■  l/<cos2kd  +  sin 


2kdl(pc)  /(pc)  ]2l 

O  S  ) 


This  factor  F  may  be  noticeably  different  from  one,  if  the  complex  specific 
acoustic  impedance  of  the  sample  material  is  not  very  close  to  that  of  the 
water  in  the  tube.  Even  though  densities  and  propagation  speeds  of  sample  and 
water  might  be  close,  as  is  the  case  for  elastomers,  a  difference  in 
attenuation  still  would  make  this  factor  unequal  to  one. 

This  proposed  correction  does  not  consider  the  effect  of  the  termination 
of  the  tube  by  the  reflector  and  the  consequent  deviation  from  the  idealized 
situation  of  an  infinite  waveguide. 
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